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PREFACE 


This  report  details  work  undertaken  under  contract  DAJA-45-87-C-005 3  to 
July  1989.  We  report  the  further  development  of  the  ungauged  forecasting 
model  MILHY.  Specifically,  new  routines  are  introduced  to  allow  discrete 
routing  in  channel  and  floodplain  zones,  and  allowance  is  made  for 
turbulent  exchange  between  channel  and  floodplain.  TKo  program  (MILHY3) 
is  applied  to  the  Fulda  watershed  in  West  Germany. 

In  addition,  a  2-dimensional  finite  element  scheme  (RMA2)  is  applied  to  a 
30  km  reach  of  the  River  Fulda.  We  conclude  that  it  should  be  perfectly 
feasible  to  couple  MILHY3  to  RMA2 :  MILHY 3  generating  the  inflow 
hydrograph  and  RMA2  predicting  the  detail  of  inundation  downstream  as 
well  as  predicting  stage  at  the  final  outflow.  This  latter  component  of 
the  study  has  been  undertaken  at  the  Hydrologic  Engineering  Center  at 
Davis,  California.  With  funds  available  at  Bristol  University  we  propose 
to  continue  this  research  and  development  work  over  the  next  18  months  to 
October  1989  under  the  aegis  of  the  contract  research  area. 

The  code  for  MILHY3,  as  developed  on  a  Sun  3/60  workstation  by  Laura 
Baird,  is  contained  in  an  appendix  to  this  report. 
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1.  IN'TRODL'CTI  jS 


Background 


This  study  relates  to  the  further  development  in  operational  model  for 
ungauged  catchment  forecasting.  The  model  used  as  the  starting  basis  for 
the  project  was  MILHY2;  a  nodel  delivered  to  Waterways  Experiment 
Station  in  1986  under  contract  DAJA-45-83-C-0029 . 


The  history  of  MILHY  development  as  an  ungauged  forecasting  model  and 
research  scheme  is  as  follows: 


MILHY :  model  for  ungauged  flow  forecasattng  using  Curve  Number  (CN) 

scheme  to  generate  runoff,  1982 

MILHY l :  adaptation  of  MILHY  under  contracts  DAJA-37-82-C-0092  and 

DAJA-37-81-C-0221  by  Dr  M  G  Anderson  to  replace  CN  scheme  by  a 
physically  based  runoff  generation  method  (finite  difference) 
1984 


MILHY2 :  development  and  validation  of  MILHY1  on  small  subcatchraent 

2 

scale  (  1  km  )  watersheds,  by  Dr  M  G  Anderson  and  Dr  S  Howes 

under  contract  DAJA-45-83-C-0029 .  Code  delivered  to  WES  in 
1986. 


MILHY 3 :  further  development  of  MILHY2  and  the  subject  of  research  In 

the  current  contract  as  specified  in  Figure  1.2  below. 

Upon  the  Initiation  of  the  current  contract  MILHY2  represented  a  fully 
working  scheme  (figures  1.1  and  1.2).  MILHY2  was  subjected  to  limited 
validation  as  shown  in  table  1.1  and  figures  1.3  and  1.4.  The  main 
conclusions  of  contract  DAJA-45-83-C-0029  relating  to  the  development  of 
MILHY2  were  that : 

(i)  the  correlation  between  predicted  and  measured  peak  discharge 
using  MILHY2  was  high  (r  -  0.91) 
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(ii)  the  time  to  peak  discharge  estimattin  was  good  using  MI1HY2  .r  = 

0.97) 

(iii)  a  comparison  of  MILKY  and  MILHY2  for  32  experimental  frames  showei 
strong  evidence  of  the  overall  improvement  achieved  by  MILHY2. 


There  was,  therefore,  a  strong  basis  for  pursuing  the  further  development 
of  MILHY2,  to  explore  the  scope  for  improvement  in  the  channel  routing 
component . 


1.2  Objectives  and  scope 

The  overall  objective  is  the  further  development  of  MILHY2  in  the 
following  specific  areas: 


(1)  employment  of  alternative  equations  for  flow  conveyance  in 
compound  channels 

(Ii)  the  improvement  In  the  handling  of  out-of-bank  roughness 

(iii)  the  validation  of  the  revised  MILKY  scheme  on  the  West  German 
watersheds  of  the  Fulda  and  Haune . 

It  has  been  possible  within  the  scope  of  the  project. to  augment  the  above 
objectives  with  that  of  evaluating  the  performance  of  a  finite  element 
model  (RMA2)  in  estimating  floodplain  Inundation.  This  work  has  been 
undertaken  at  the  Hydrologic  Engineering  Center,  Davis  California,  and 
will  continue,  with  available  funds  at  Bristol  University,  for  a  further 
eighteen  months  to  October  1990. 


The  above  objectives  can  be  directly  translated  to  research  questions: 
(1)  can  MILHY3  accommodate  improved  channel  routing  procedures? 


(Ii)  can  MILHY3,  as  presented  schematically  In  figure  1.5  satisfy 
the  restricted  data  needs  for  an  ungauged  model? 
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Table  1.1: 


Comparison  of  catchment  characteristics  which  are 
require-!  by  the  unit  hydrograph  procedure 


Area 

(km2 

Difference 
in  elevation 
(ra) 

length  of 

Main  channel 
(ka) 

W-2 

North  Danville 
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0.6 

79.3 

1  #  > 
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Figure  1.4  :  Comparison  of  percentage  time  to 

peak  discharge  error  of  MILHY2  and 
HYMO,  32  experimental 
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RESEARCH  DESIGN' 


2.  i  Introd action 

The  overall  research  design  of  the  M1LHY  project  as  developed  bv  Or. 
Anderson's  group  at  Bristol  University  over  the  iasc  six  years  is  shown 
in  figure  2.1.  The  initial  decision  regarding  MILHY1  related  to 
utilisation  of  finite  difference  methods  for  runoff  generation. 

Subsequent  research  identification  suggested  the  need  to  examine 
alternatives  for  compound  channel  modelling  and  it  is  this  development 
that  is  the  major  constituent  of  MILHY3.  However,  more  general  issues 
are  raised  here  in  the  context  of  the  interaction  of  hydraulic  and 
hydrologic  schemes,  and  their  respective  suitability  for  ungauged 
inundation  modelling. 

With  the  development  of  MILHY3  to  larger  watershed  scales  (  1000  km") 

the  potential  for  hydraulic  handling  of  this  problem  is  deemed  worthy  of 
Investigation.  As  we  will  discuss  In  Chapter  7  of  this  report,  little 
work  has  been  done  on  the  application  of  2-dimensional  finite  element 
models  for  river  applications  at  this  scale,  let  alone  exploring  the 
suggestion  we  make  in  our  overall  research  strategy  for  examination  of 
the  possibility  of  coupling  hydrological  (MILHY3)  and  hydraulic  schemes 
(RMA2)  -  see  figure  2.1. 

2 . 2  Research  design  for  the  current  Investigation 

Any  model  design  Is  essentially  a  two  dimensional  matrix  of  components. 
This  la  Illustrated  In  general  terras  in  figure  2.2.  Decisions  have  to  be 
made  In  the  two  principal  areas:  (l)submodel  Inclusion,  and  (2) 
resolution  of  the  selected  submodel.  The  area  of  submodel  Inclusion  in 
general  terms  has  been  left  unaltered  from  earlier  versions  of  MILHY.  It 
Is  the  latter  decision  area  that  has  proved  the  focus  of  the  current 
research.  In  particular,  we  have  sought  to  examine  the  effects  of 
changes  in  the  channel  routing  submodel,  In  the  context  of  examining  and 
implementing  alternative  models  for  compound  channel  flow  conditions.  In 
addition,  a  somewhat  lesser  effort  has  been  expended  In  an  examination  of 
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Che  precipitation.  Figure  2.2  shows  the  scope  of  che  submodel  resolution 
development.  An  important  consideration  here  relates  to  achieving  a 
submodel  resolution  that  is  considered,  broadly  at  least,  to  be 
consistent  between  all  submodels.  This  is  an  issue  to  which  we  will 
return;  submodel  resolution  cannot  of  course  be  divorced  from  the 
resolution  of  the  user  supplied  information  (figure  1.1).  Varying, 
perhaps  unavoidable,  resolution  in  the  user  supplied  information  may  be 
considered,  potentially  at  least,  to  have  significant  ramifications  for 
submodel  and  model  overall  performance.  Thus  user  supplied  data 
resolution  cannot,  and  should  not,  be  divorced  from  model  formulation, 
design  and  validation.  Regrettably,  in  many  cases,  this  association  is 
not  made.  As  we  will  discuss  in  later  sections,  this  concept  is  central 
to  the  issue  of  design,  implementation  and  validation  of  MILHY3. 
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Figure  2.1  :  Bristol  University  overall  research  design 
for  the  MILHY  project 


Figure  2.2  :  Component  logic  structure  for  MILHY3 
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1  3.  IDENTIFICATION  OF  KEY  PARAMETERS  AND  PROCESSES  -  )R  C AN 


To  achieve  the  main  objective  of  improving  the  predictive  capabil itv  of 
the  conveyance  section  of  MILHY,  it  was  important  firstly  to  identify  the 
most  important  processes  and  parameters.  An  optimal  method  of 
identification  wouLd  be  to  undertake  a  sensitivity  analysis  of  an 
existing  model  of  two-stage  flow. 

3 .  t _ Difficulties  of  modelling  two-stage  channels 

Two-stage  channels  consist  of  a  main  channel  and  adjoining  floodplains 
which  are  subject  to  inundation.  Water  on  the  floodplains  may  oe  either 
stationary,  where  the  floodplains  act  as  stores  of  water,  or  flowing  when 
the  floodplains  act  as  another  channel. 

3.1.1  The  complexity  of  physical  processes 

Two-stage  channels  are  a  complex  three-dimensional  system.  The  inclusion 
of  the  floodplain  system  is  not  simply  a  matter  of  extending  the 
cross-sectional  area  of  the  main  channel.  As  Bhowmik  and  Demissie  (1982) 
have  shown,  the  carrying  capacity  of  a  two-stage  channel  is  not  directly 
proportional  to  cross-sectional  area.  Figure  3.1  illustrates  the 
theoretical  line  of  proportionality  between  area  and  discharge,  and  the 
relationship  developed  from  field  data  collected  from  five  rivers  in  the 
USA.  From  figure  3.1  it  is  possible  to  conclude  that  two-stage  channels 
may  not  be  considered  as  a  single  system. 

Bhowmik  and  Demissie's  data  also  confirmed  flume  studies  by  Rajaratnam 
and  Ahmadi  (1979),  that  there  is  interaction  between  the  water  on  the 
floodplain  and  water  in  the  main  channel.  Rajaratum  and  Ahmadi  (1979) 
found  discontinuity  in  the  velocity  fields  of  the  main  channel  and  bed 
shear  at  the  boundary  between  the  channel  and  floodplain.  Figure  3.2 
shows  the  stage/ velocity  relationship  in  Salt  Creek,  USA,  and  illustrates 
the  discontinuity  of  velocity  in  the  main  channel. 


A. 


T 


■! 

. 

i 

i 


P 


15 


Discharge  in  floodplain  /  total  discharge 


Figure  3.1  :  Relationship  between  the  carrying  capacity 
and  floodplain  area  (after  3howmik  and 
Demissie,  1982) 


A. 


SALT  CREEK 


V  Cfps) 

Figure  3.2  :  Stage/ velocity  relationship  in  a 

two-stage  channel  (after  Showmik  and 
Demissie,  1982) 
a  -  floodplain 
b  -  combined 
c  -  channel 
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The  reduction  in  the  velocity  or  flow  in  Che  tain  channel  is  rearmin'--, 
conditions  are  reached  suggested  that  there  was  transverse  mass  transfer 
between  the  floodplain  and  the  main  channel,  effectively  retarding  flows 
in  the  main  channel  and  accelerating  them  in  the  floodplain.  This 
momentum  transfer  may  be  envisaged  to  occur  through  the  action  if 
turbulent  shear  stresses,  first  recognised  and  photographed  by  Sellin 
(1964). 

From  figure  3.2,  it  can  be  seen  that  Che  minimum  channel  velocity  occurs 
when  the  floodplain  stage  elevation  is  approximately  35%  of  Che  main 
channel  stage,  and  that  as  floodplain  stage  increases  the  system 
converges  to  the  composite  value.  Bhowmik  and  Demlssie's  results 
suggest,  therefore,  that  the  behaviour  of  the  two-stage  channel  depends 
on  the  depth  of  Inundation  of  the  floodplain. 

Plan  geometry  must  also  be  included  when  attempting  to  expose  the 
processes  accive  in  two-stage  channels.  Varying  cross-sectional 
geometries  associated  with  meandering  main  channels  cause  floodplain 
inundation  at  different  points  along  a  reach.  Wolman  and  Leopold  11957) 
showed  chat  the  return  period  of  floodplain  inundation  increased  as  the 
width/depth  ratio  of  the  main  channel  increased.  The  downstream  and 
orthogonal  floodplain  slopes  then  determine  if  the  floodplain  water 
rejoins  the  main  channel  or,  as  Fread  (1976)  suggests,  routes  downstream 
along  a  different  path  to  the  main  channel. 

Toebes  and  Sooky  (1967)  in  a  series  of  flume  experiments  showed  that  the 
momentum  transfer  between  floodplain  and  main  channel  flows  are 
exacerbated  where  the  floodplain  flow  Is  at  an  angle  to  the  main  channel 
flow,  primarily  In  meandering  flows.  Chang  (1983)  showed  that  the 
Increased  energy  expenditure  in  a  mature  meandering  channel  Is  due  to: 

1)  Internal  fluid  friction  caused  by  transverse  circulation 
(secondary  currents) 

li)  boundary  resistance  associated  with  transverse  shear 
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With  data  from  Poudre  Supply  Canal,  Fort  Collins,  Colorado,  Chang  was 
able  to  conclude  that  where  the  depth/radius  is  high  or  the  roughness 
is  low,  the  transverse  circulation  energy  losses  may  be  greater  than 
chose  associated  with  the  primary  or  longitudinal  flow. 

3.1.2  Modelling  alternatives 


Research  into  modelling  of  two-stage  channels  has  included  both  hvdraul 
and  hydrologic  approaches.  Hydraulic  engineers  solve  both  the 
conservation  of  mass  and  a  simplified  form  of  the  conservation  of 
momentum  (equations  3.1  and  3.2),  whilst  hydrologists  use  only  the 
conservation  of  mass  (equation  3.3)  and  a  relationship  between  flow  and 
storage . 


Conservation  of  mass: 

2(_AV)  +  9a  -  0 

t 


3.1 


Conservation  of  momentum- 

"dv  +  vSv  +  g  (3h  +  S  )  -  0 

57  3x  (Jx  0 

Conservation  of  mass: 

I  -  0  -  As/At 


3.2 


3.3 


A  *  area 

V  »  velocity 

x  -  longitudinal  axis 
g  -  gravitational  acceleration 

h  -  stage 


S  »  friction  slope 

o 

t  -  time 

I  •  Inflow 

0  •  outflow 

S  »  storage 
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Hydrologic  approaches  have  been  limited  to  one- dimension ,  whilst 
hydraulic  models  include  both  one  and  two  dimensions  and  even  a  protofoe 
three-dimensional  turbulence  model  (Krisnappen  and  Lau,  1966). 

One-dimensional  approaches  (solving  the  St.  Venant  equations  of  flow  in 
hydrauLic  models)  either- 

1)  treat  the  channel/f Loodpiain  cross-section  as  a  single  system 
and  average  boundary  roughness  and  velocities 

2)  or  treat  the  floodplain  as  an  area  for  storing  water  only,  for 
example,  the  Hydrologic  Engineering  Centre  model  (HEC-1,  1981) 

3)  or  divide  the  cross-section  into  homogeneous  segments  of  rlow 
but  do  not  consider  momentum  transfers  between  these  segments, 
for  example  Tingsanchali  and  Ackermann  (1976). 

The  two-dimensional  process  of  momentum  transfer  between  the  main  channel 
and  floodplain  segments  can  only  be  modelled  using  a  two-dimensional 
hydraulic  approach  using  the  Reynolds  equations.  However,  empirical 
relationships  have  been  derived  for  the  relationship  between  the  boundary 
shear  stresses,  (apparent  stresses  produced  by  momentum  transfer)  and  the 
cross-sectional  geometry  of  the  two-stage  channel  (Knight  and  Demetriou, 
1983).  Pasche  and  Rouve  (1985)  have  developed  a  one-dimensional  model 
Incorporating  these  boundary  shear  stresses  and  a  hydraulically-based 
velocity  distribution.  Incorporation  of  these  boundary  shear  stresses  is 
more  fully  investigated  in  section  4. 

Modelling  of  the  process  of  momentum  transfer  has  been  incorporated  in 
two-dimensional  hydraulically-based  finite-difference  and  finite-element 
models.  They  usually  employ  one  of  the  three  methods  below  to  quantify 
the  effects  of  momentum  transfer: 

1)  compute  the  force  to  provide  equilibrium  in  each  segment  of 
flow  (apparent  boundary  shear  force) 
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2)  compute  effective  friccion  factors  for  each  segment 

3)  compute  the  iteration  between  a  shear  layer  and  the  velocity 
profile  (turbulence  model) 

I  An  investigation  of  these  models  and  the  application  of  a  finite-element 

model  to  a  two-stage  channel  is  reported  in  Section  7. 


3.1.3  MILHY  :  present  frictional  capability 


The  handling  of  friction  has  remained  unaltered  in  all  versions  of  MILHY 
where  it  is  handled  as  the  Manning  'n'  coefficient,  utilized  in  the 
Manning  equation  (equation  3.4)  to  compute  the  stage/discharge 
capability: 


i 

l 


V 


R1/3  .  S1 2 3 4 5 6 7 8 * 10^ 
n 


3.4 


The  cross  section  is  divided  up  by  the  user  such  that  each  segment  of  the 
section  is  fractionally  homogeneous  and  hence  one  value  of  'n'  per 
segment  Is  entered  to  compute  the  stage/discharge  relationship.  In 
selecting  the  most  appropriate  'n'  value  for  each  segment,  the  user  must 
consider  the  following  factors 

1)  surface  roughness 

2)  vegetation 

3)  channel  Irregularity 

4)  channel  alignment 

5)  silting  and  scouring 

6)  obstructions 

7)  size  and  3hape  of  channel 

8)  stage  and  discharge 

9  seasonal  change 

10)  suspended  material  and  bed  load 


i 

I 
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Xc  is  a  difficult  task  to  select  one  empirical  'o'  value  chat 
successfully  describes  a  particular  channel.  Manning's  'n'  attempts  to 
incorporate  many  of  the  physical  processes  active  in  the  river  channel 
system  as  well  as  the  effects  of  boundary  friction.  The  aim  of  the 
sensitivity  analysis  reported  below  is  to  Identify  which  of  the  physical 
processes  active  in  a  two-stage  channel  are  dominant.  In  later  sections, 
alternative  methods  to  the  Manning's  'n'  coefficient  of  handling  these 
processes  are  investigated. 

At  present,  the  only  additional  complexity  to  he  Manning's  'n'  handling 
of  friction  incorporated  in  MILHY,  is  an  algorithm  (equation  3.5)  to 
reduce  the  'n'  value  with  increasing  stage. 

n'  =  n  -  0.0025R  3.5 

If  the  dominant  process  active  in  a  channel  is  boundary  roughness  then 
this  algorithm  will  improve  the  prediction  of  the  carrying  capacity  of 
the  cross-section.  In  the  main  channel  as  stage  increases  then  the  area 
of  fLow  will  increase  faster  than  the  wetted  perimeter,  thus  reducing  the 
retarding  effects  of  friction  along  the  wetted  boundary  (SCS,  1954).  On 
the  floodplains  too,  Manning's  'n'  may  decrease  as  the  depth  of 
Inundation  increases  and  the  frictional  effects  of  vegetation  become  less 
important.  Table  3.1,  taken  from  Chow  (1959),  illustrates  this  decline 
for  pasture  and  meadows,  typical  in  the  floodplains  of  the  River  Fulda 
catchment  in  Wer  •  Germany. 

Petryk  and  Bosmajlan  (1975)  showed,  however,  that  this  Is  an 
over-simplification  of  the  frictional  effects  of  vegetation.  When  the 
vegetation  is  submerged  then  'n'  will  decrease  with  increasing  stage,  but 
below  the  top  of  the  vegetation  then  there  is  a  complex  relationship 
between  the  vegetation  density  and  Manning's  'n'.  Petryk  and  Bosmajlan 
(1975)  used  equation  3.6  to  calculate  the  change  in  'n'  with  depth. 
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3.6 


TabLe  3.1 


Mannings  ' n'  value  for  pasture  and  meadow  floodplains 
from  Chow  ( 1939) 


Depth  of  Inundation  Mannings  'a'  value 


( feet) 

Pasture 

Meadow 

Less  than  l 

0.05 

).[■: 

1-2 

0.05 

0.01 

2-3 

0.04 

O 

o 

3-4 

0.04 

0.0< 

Over  4 


0.04 


0.0 


where  n  =  Manning  'n'  value  with  vegetation  effects 

D 

=  vegetation  drag  coefficient  dependent  on  vegetation  type 
L  =  length  of  reach 

.V  =  projected  area  of  the  ith  plant 
A  =  cross-sectional  area  of  flow 

petryk  and  Bosmajian  (1975)  scheme  would  not  be  suitable  for  the  ungauged 
application,  however  it  does  show  the  inability  of  one  Manning's  'n' 
value  to  represent  the  effects  of  all  the  conditions  listed  above.  In 
two-stage  channels  this  is  even  more  of  a  problem  as  the  dominant 
processes  change  as  the  floodwave  overtops  the  main  channel  and  floods 
out  over  the  floodplains. 

The  n  reduction  algorithm,  equation  3.5,  assumes  that  boundary  roughness 
is  the  dominant  process;  in  the  two-stage  channel  this  is  not  the  case. 
As  stage  increases  enough  so  that  the  floodplain  is  inundated,  not  only 
is  there  a  rapid  increase  in  the  wetted  perimeter  but  turbulent  eddies 
between  the  main  channel  and  floodplains  occur.  The  frictional  effects 
of  these  eddies  are  much  greater  than  the  increase  in  boundary  friction 
(Pasche  and  Rouve ,  1985). 

However,  in  application  of  MILHY2  to  two-stage  channels  in  the  River 
Fulda  catchment,  a  more  immediate  problem  occurred.  As  floodplain 
inundation  depths  increased,  the  hydraulic  radius  Increased  such  that  n' 
in  equation  3.5  became  negative.  An  example  of  this  problem  is  found  in 
Table  3.2,  illustrating  a  typical  stage/discharge  relationship  produced 
by  MILHY2  under  out-of-bank  conditions.  For  this  reason  the  algorithm 
has  been  removed  from  MILHY3  and  is  not  included  in  any  of  the 
simulations  reported  here. 

The  atm  of  the  work  reported  here  is  therefore  to  retain  the  Manning  'n' 
coefficient  to  Incorporate  the  boundary  friction  effects  on  the 
floodplain  and  all  of  the  ten  factors  identified  earlier  for  the  ln-bank 
bank  channel.  Incorporation  of  the  additional  frictional  effects  in  the 
two-stage  are  Investigated  in  section  4  and  5. 
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Table  3.2 


MILHY2  ROUGHNESS 
REDUCTION  FORMULA: 
a’  =  n  -  0.0025R 


Rating  curve  valley 

section 

Water 

Flow 

Flow 

surface 

area 

rate 

CFS  (x 1 :  ) 

Elev . 

sq,  ft. 

591.38 

86.4 

0 

595.45 

299.6 

1 

599.02 

573.3 

3 

602.59 

1030.0 

7 

606.16 

4106.6 

12 

609.73 

8436.7 

63 

613.30 

12920.8 

256 

616.86 

18078.9 

356  1  5 

620.43 

23984.5 

-54 

624.00 

30582.6 

-327 
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3 . 2 _ Selection  of  a  two-stage  conveyance  model 

From  section  3.1  it  would  seem  that  a  successful  conveyance  model  must 
incorporate : 

1)  plan  and  cross-sectional  geometries 

2)  momentum  transfer  of  flow  between  floodplain  and 
main  channel  segments 

An  ideal  solution  therefore  would  be  a  two-dimensional  hydraulic  model. 
However,  as  Section  7  will  show,  such  models  have  not  been  applied  to  the 
scale  of  the  reach  under  investigation  here  (i.e.  greater  chan  10  km  in 
length).  A  one-dimensional  approach  was  therefore  accepted. 

From  the  literature  it  seemed  that  most  of  the  investigations  into 
two-stage  channel  conveyance  have  concentrated  on  the  transfer  of 
momentum  across  the  section  and  techniques  for  its  inclusion  are  well 
established.  It  was  seen  as  a  priority,  therefore,  to  select  a 
one-dimensional  model  which  would  permit  concentration  of  the  sensitivity 
analysis  on  plan  geometry  processes  and  parameters.  Hence,  the  model 
selected  was  a  state-of-the-art  model  developed  by  Ervine  and  Ellis 
(1937),  which  considers  the  three-dimensional  system  using  a  series  of 
analytical  equations. 

3.2.1  The  Ervine  and  Ellis  Model 


Ervine  and  Ellis'  model  allows  a  meandering  plan  geometry  to  be  modelled 
by  dividing  flow  into  three  segments,  shown  on  figures  3.3a  and  b,  and 
defined  as: 

1)  main  channel  flow 

2)  floodplain  flow  contained  within  the  meander  belt  of  the 
main  channel 


3) 


floodplain  flow  outside  of  the  meander  belt 
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For  each  segment,  the  energy  loss  is  computed  and  hence  the  mean  velocity 
for  each  segment  and  discharge  cotal  are  calculated.  Ervine  and  Ellis 
firstly  identified  the  main  sources  of  energy  loss  of  each  segment  of 
flow  and  then  brought  together  a  series  of  geometric  and  fictional 
relationships  to  explain  them. 

Main  channel  energy  losses 


Ervine  and  Ellis  considered  there  to  be  four  possible  sources  of  energy 
loss 

a)  frictional  losses  at  the  boundaries 

b)  transverse  currents  (secondary  currents)  at  meander  bends 

c)  turbulent  shear  stress  (momentum  transfer  to  floodplains) 

d)  pool/riffle  sequences  causing  head  loss  at  low  flows 

Ervine  and  Ellis  chose  to  omit  the  turbulent  shear  stresses  and 
pool/riffle  losses  in  their  computation.  Shear  stresses  were  omitted 
because  three-dimensional  interpretations  of  established  techniques  (e.g. 
Knight  and  Demetriou  (1983))  are  still  under  investigation  by  Willetts 
(see  Ervine  and  Ellis  (1987)).  Pool/riffle  losses  are  considered  less 
important  in  times  of  overbank  flow,  when  bed  form  effects  are  usually 
flooded  out. 

Floodplain  energy  losses 

Two  sources  of  loss  were  identified: 

a)  frictional  losses  at  the  boundaries 

b)  expansion  and  contraction  losses,  shown  in  figure  3.4,  where 
flow  orthogonal  to  the  main  channel,  suddenly  expands  as  it 
drops  into  the  channel,  and  contracts  as  it  re-enters  the 
floodplain  region 
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EXPANSION! 


C0MTRAC.7I0W 


ii 


Fig .  3.4 


3.2.2  Quantifying  the  energy  losses 
Main  channel  energy  losses 


a)  Friction 

Head  losses  due  to  friction  are  computed  over  a  meander 
wavelength  (  r^m)  as : 


where  f^  is  the  Colebrook-White  friction  factor,  given  by 


b)  Transverse  currents 

Head  losses  due  to  secondary  currents  at  meander  bends  are 
computed  using  a  simplified  method  developed  by  Chang  (1983). 
Chang  used  a  mean  transverse  current  'secondary  current)  velocity 
because  over  a  meander  amplitude,  the  velocity  varies  from  a 
maximum  at  the  apex  to  a  theoretical  zero  at  the  cross-over 
thalweg.  Chang  Ignored  the  effects  of  superelevation,  where 
centrifugal  forces  cause  the  water  level  on  the  outside  of  the 
bend  to  be  higher  than  those  on  the  Inside;  (Yen,  1967,  showed 
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chat  in  two-stage  flow,  superelevation  effects  are  suppressed  by 
the  head  of  water  above  the  -nain  channel). 


Head  loss  due  to  transverse  currents  computed  aver  a  meander 


wavelength  is 

given  by: 

hi  -* 

Vfe.  T 

1-0*  fc] 

— 

IT 

u  O-SfeS 

— 1 

V 

UJ  la,  ) 

Floodpl 

ain  within  the 

meander  belt 

width 

Friction 
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As  in  the  main  channel,  the  total  frictional  head  loss  along  a 
meander  wavelength  is  described  by: 

j  WmVm-rVmBc. 

where  the  last  term  is  the  wetted  area. 

b)  Expansion  losses 

Assuming;  y^  -  +•  h  3.11 

the  head  loss  due  to  expansion  of  floodplain  flow  into  the  main 
channel  over  a  meander  wavelength  is  given  by: 

he  -  rWl-yfV  V 

^  ye/  2.3 

where  Q  is  the  average  mean  angle  of  the 
main  channel  over  a  meander  wavelength. 
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3.12 


floodplain  flow  to  the 


A. 
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c)  Contraction  losses 

The  total  head  loss  due  to  contraction  if  floodplain  flow  leaving 
the  main  channel  (illustrated  in  fig.  3.4)  is  given  by 

Hc*  CL.  .  si n.2  (  rVrrO 

^9 

where  is  a  loss  coefficient,  generated  by  Yen  and  Yen  (1984), 
and  is  a  function  of: 

i)  the  density,  specific  weight  and  kinematic  viscosity 
of  the  flow 

ii)  meander  wave  length  and  amplitude,  the  mean  angle  of 
incidence  of  floodplain  flow  to  the  main  channel,  the 
valley  width,  valley  slope,  floodplain  roughness,  and  the 
width  and  depth  of  the  main  channel 

iii)  discharge  and  shape 

Yen  and  Yen  (1984)  using  data  collected  from  flume  experiments 
computed  Che  total  loss  coefficient  after  flow  had  been  subjected 
to  expansion  and  contraction.  Then  assuming 

where  C  3  total  loss  coefficient  3.14 

C£  -  loss  coefficient  due  to 
expansion 

»  loss  coefficient  due  to 
contraction 


c  »  4  +  Cl 


and  C£ 


2 


as  in  3.12 


3.15 
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the  contraction  loss  coefficients  were  compute!,  is  shown  in 

Table  3,3. 

Yen  and  Yen  consider  that  these  coefficients  should  be  created  as 
upper  limits  because  the  channeL  sidewalls  in  their  flume 
experiments  were  vertical  rather  than  a  more  realistic  trapezoidal 
slope. 

Floodplain  flow  outside  the  meander  belt 
a)  Friction 


In  the  floodplain  outside  the  meander  belt,  flow  is  considered  to  be 
uniform,  therefore  the  friction  slope  is  given  by 


Combining  all  the  head  loss  equations,  Ervine  and  Ellis  (1987) 
obtained 


i)  for  the  main  channel 


(t sV'-Wjv'/  I 

Y#Vfa+  2.0*fcl 

UA  r  Jt<j  [ 

.  0- 5fc5  +  Vfc  J 

(j )  • 

■  3  Xfm 

*9 

ii)  for  the  floodplain  Inside  the  meander  belt  width 


£ 


* 


Table  3.3  :  Contraction  Loss  Coefficients 
(after  'fen  and  Yen,  i984) 


0.0 

0.1 

0.2 

0.3 

0.4 

0.5 

0.6 

0.7 

0.8 

3.9 

0.3 

0.48 

3.45 

0.41 

0.36 

0.29 

9.21 

9.13 

0.07 

9.0 
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iii)  for  the  floodplain  outside  the 


ier  belt  .-id: 


3.  '.9 


alongside  with  total  discharge,  using  the  continuity  equation  which  is 
given  by 


Q  =  Vc  (Bch) 


Vfl  (yf  V  +  Vf2  yf  (Wt  -  V 


3 . 3  Sensitivity  Analysis  of  the  Ervine  and  Ellis  Model 


3.:n 


The  objective  of  undertaking  a  sensitivity  analysis  of  Ervine  and  Ellis' 
model  was  to  identify  the  physical  processes  controlling  the  velocity  and 
discharge  predictions.  Once  identified  the  most  appropriate  method  of 
incorporating  these  processes  into  the  MILHY  scheme  can  be  investigated. 

3.3.1  Sensitivity  analysis  design 

Analysis  of  equations  3.7  to  3.16  and  3.20  identified  five  groups  of 
parameters  which  controlled  the  processes  identified  and  modelled  by  Ervine 
and  Ellis.  These  five  groups  are:- 

i)  slope 

ii)  plan  geometry  (channel  width,  floodplain,  meander  belt  width  and 
radius  of  curvature 

iii)  depth  of  flow  (in  channel  and  floodplain  segments) 

iv)  sinuosity  (sinuosity  and  angle  of  Incidence  of  floodplain  flow 
to  the  main  channel) 

v)  friction  (for  consistency  with  later  MILHY  analysis  the  sensitivity 
of  the  model  to  Manning  'n'  was  used,  utilizing  the  conversion 
equation  3.21,  below): 

f  »  8gn~ 

R 1  /  3 


3.21 
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Each  of  these  five  groups  were  investigated  individuai.lv  bv  varying  eath 
group  independently  by  a  systematic  30’;  and  5"  reduction  and  5’;  and  30’; 
increase  in  parameter  values. 

The  model  was  applied  to  a  hypothetical  reach  because  the  sinuosity 
variables  only  allow  for  a  constant  amplitude  and  wavelength  values  for  the 
main  channel  throughout  the  reach  length.  However,  all  parameter  values 
were  taken  from  topographic  information  collected  from  the  River  Fulda,  West 
Germany.  The  values  applied  are  shown  in  Table  3.4  and  the  velocity 
predictions  for  each  flow  segmens  and  discharge  predictions  by  the  model  for 
these  values  are  also  shown.  Observed  stage/discharge  relationships  from 
the  River  Fulda  confirm  that  these  model  predictions  give  realistic 
discharge  values. 

3.3.2  Results 

The  results  from  the  sensitivity  analysis  are  tabulated  in  Tables  3.3  to 
3.9,  and  show  the  percentage  deviation  from  .ae  computed  values  tabulated  in 
Table  3.4.  Below  is  an  analysis  of  the  velocity  predictions  by  considering 
each  of  the  sources  of  head  loss  identified  by  Ervlne  and  Ellis: 

t)  Frictional  losses  are  modelled  in  all  three  flow  segments  and  Table 
3.5  to  3.7  show  that  variation  in  the  frictional  parameter  values 
cause  the  largest  variation  in  the  velocity  predictions.  However,  in 
the  main  channel  the  Darcy-We Isbach  friction  factor  is  also  linked  to 
the  modelling  of  the  transverse  (secondary)  circulation.  From  the 
first  term  in  equation  3.9,  it  can  be  seen  that  as  the  friction 
factor  decreases,  head  losses  from  the  transverse  currents  decrease, 
and  when  the  friction  factor  increases  head  Losses  are  increased. 
Therefore  the  velocity  variations  shown  in  Table  3.5  incorporate  both 
friction  head  losses  and  transverse  circulation  losses. 

li)  The  transverse  circulation  in  the  main  channel  can  be  attributed  to 

friction  (as  noted  above)  and  the  ratio  of  hydraulic  radius  to  radius 
of  curvature.  This  ratio  is  included  in  the  geometry  variation 
reported  on  Table  3.5,  which  shows  that  the  velocity  predictions  are 
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Table  3.4:  Parameter  specification  for  hypothetical  reach 


SI  units 

3ed  slope  0  .  3007 

Sinuosity  1  .  3 

Hydraulic  radius  2.5 

Radius  of  curvature  125.0 

Width  of  meander  belt  175.0 

Total  floodplain  width  300.0 

Channel  width  30.3 

Friction  channel  (f)  0.371 

Friction  floodplain  1  0.356 

Friction  floodplain  2  0.356 

Depth  channel  3.5 

Depth  flood  plain  0.5 

Angle  of  flood  plain  flow  to 

channel  (radians)  0.785 

Contraction  loss  coefficient  0.47 


Results 


Main  channel  velocity  1.205 
Floodplain,  area  1  velocity  0.360 
Floodplain,  area  2  velocity  0.278 
Discharge  157.2 
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Table  5.5:  Channel  Velocity  Results 
(”  deviation  from  origin  velocity) 


”  Change  in 
variable 

Decrease 

30s 

Decrease 

5% 

Increase 

5% 

Increase 

3  V; 

Slope 

-19 

+  3 

+  13 

Channel 

friction 

+•50 

+  5 

-  4 

-21 

Geometry 

-  5 

LTi 

O 

1 

+  i 

+  3 

Sinuosity 

+20 

+3 

-a 

-  12 
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Table 


3.6:  Flood  plain  (area  1)  velocity 

(’(  deviation  from  origin  velocity) 


results 


X  Change  in 
variable 

Decrease 

30"< 

Decrease 

5? 

Increase 

5  X 

Increase 
3  v- 

Slope 

-19 

-  2 

+  3 

+  13 

Friction 

+13 

+  5 

-  4 

- 7  7 

Geometry 

-  4 

-  0.5 

+  1 

+  2 

Sinuosity 

+  i 

0.0 

0.0 

-  1 

Contraction 
coef  f icient 

0.0 

0.0 

0.0 

0.0 
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Table  3.7:  Floodplain  (area  2)  velocity  results 
(Z  deviation  from  origin  velocity) 


X  Change  in 
variable 

Decrease 

30? 

Decrease 

5 \ 

Increase 

5'< 

Increase 

3:; 

Slope 

-19 

_  +> 

+  3 

+  13 

Floodplain 

friction 

+25 

+  5 

-  5 

-23 

Table  3.8:  Flow  depth  effects  on  velocity  and  dischsrg 
(*  deviation  from  origin) 


Channel  Flood  Plain  Dischar 

Velocity  Velocity 

Depth  Area  1  Area  3 

( metres) 


0.33 

2.31 

-12 

0.475 

3.325 

-  2 

0.525 

3.675 

+  4 

0.665 

4.635 

+15 

0.33 

3.5 

- 

0.475 

3.5 

- 

0.525 

2.5 

- 

-33 

-23 

-29 

-20 

-20 

-19 

-15 

-16 

-16 

-  5 

-  5 

-  4 

-13 

-19 

-14 

-  2 

-  2 

__  1 

+  2 

+  3 

+-  2 

+  15 

+  15 

+ 1 6 

0.665 


3.5 
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Table  3.9:  Discharge  results 
(%  deviation  from  origin  discharge) 


"  Change  in 
variable 

Decrease 

30% 

Decrease 

5% 

Increase 

5” 

Increase 

V’.: 

Slope 

-19 

-  2 

+  3 

+  1  3 

Channel 

friction 

+35 

+  4 

-  3 

-17 

Floodplain 

friction 

+  15 

+  2 

-  1 

-  9 

Geometry 

-28 

-  4 

+  4 

+27 

+  14  +  2 


Sinuosity 


8 


9 
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not  sensitive  to  geometric  variation  in  the  channel.  As  noted  ibove, 
however,  the  model  is  sensitive  to  the  frictional  aspects  of  the 
transverse  circulation. 

iii)  Sinuosity  changes  generate  significant  variability  in  the  channel 

velocity  results  (Table  3.5).  From  equation  3.17  it  can  be  seen  that 
the  sinuosity  term  is  used  to  calculate  channel  length  in  both  the 
frictional  head  loss  and  transverse  circulation  computations.  For 
the  main  channel,  therefore,  the  model  can  be  interpreted  as  being 
sensitive  to  channel  length. 

On  the  floodplain  within  the  meander  width  belt  (area  1),  Table  3.6 
shows  the  velocity  predictions  are  not  sensitive  to  sinuosity 
variations.  From  equation  3.18  it  can  be  seen  that  sinuosity  is 
utilized  to  compute  the  flow  path  length  and  the  angle  of  incidence 
of  floodplain  flow  to  main  channel  flow  used  in  the  calculation  of 
the  expansion  and  contraction  head  loss.  From  Table  3.6  it  would 
seem  reasonable  to  conclude  that  because  of  the  linear  flow  path  the 
velocity  predictions  are  not  affected  by  the  length  of  the  path,  and 
it  is  not  necessary  to  include  the  angle  of  incidence  of  floodplain 
flow  in  the  modelling  of  expansion  and  contraction  head  losses. 

Table  3.6  also  shows  that  the  exact  value  of  the  Yen  contraction  loss 
coefficient  need  not  be  of  concern  to  the  modeller. 

iv)  The  effect  of  slope  variations  on  velocity  predictions  were  only 

significant  where  variation  was  large  (+  30%),  seen  in  Tables  3.5  to 
3.7.  For  the  purposes  of  this  sensitivity  analysis  the  frictional 
slope  (SQ)  was  assumed  to  be  parallel  to  the  bed  slope,  S,  hence 
uniform  flow  was  assumed.  However,  the  sensitivity  analysis  tested 
the  effects  of  frictional  slope  variation  on  the  velocity  and 
discharge  predictions.  The  effects  of  different  slopes  in  the  main 
channel  and  floodplain  areas  are  not  directly  included  in  the  Ervine 
and  Ellis  model. 

v)  The  Impact  of  variation  in  the  depth  of  flow  on  the  floodplain 

velocity  results  are  shown  in  Table  3.8.  Equation  3.17  shows  that 
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the  main  channel  depth  is  incorporated  in  che  velocity  computation  so 
hydraulic  radius,  and  analysis  of  intermediate  computations  in  che 
analysis  shows  it  is  the  hydraulic  radius  in  the  frictional  head  Voss 
computation  to  which  the  velocity  results  are  sensitive. 

On  the  floodplain  within  the  meander  belt,  equation  3.19  shows  it  is 
the  ratio  of  floodplain  to  channel  depth  that  is  utilized  to  compute 
expansion  and  contraction  head  losses  (see  the  lower  four  results  on 
Table  3.8).  However,  the  velocity  predictions  for  floodplain  area  1 
are  identical  to  those  in  floodplain  area  1  and  as  the  headloss  in 
area  2  is  entirely  attributable  to  friction  (equation  3.19),  it  would 
seem  that  the  velocity  variations  in  areal  are  due  to  che  same 
frictional  effects,  and  not  due  to  expansion  and  contraction  losses. 

Tables  3.8  and  3.9  iliustrace  the  effects  of  variability  in  che  five  groups 
identified  on  the  discharge  predictions  computed  using  equation  3.20. 
Geomecry  is  the  only  group  to  create  additional  influence  on  the  discharge 
predictions,  over  those  on  che  velocity  resuLts  reported  above.  The 
geometry  variables  effectively  weight  the  velocity  results  for  each  flow 
area  based  on  their  cross-sectional  area,  to  give  total  discharge. 

3.3.3  Conclusions 


From  the  analysis  of  the  results  above,  it  is  possible  to  make  several 
conclusions 

i)  The  Ervine  and  Ellis  scheme  is  highly  sensitive  to  che  Darcy-Weisbach 
friction  factor. 

ii)  The  model  is  sensitive  to  the  depth  of  inundation  (incorporated  in 
the  computation  of  frictional  headlosses)  in  all  flow  areas. 

lit)  The  sinuosity  of  the  main  channel  is  Important  in  determining  the 
length  of  the  flow  path  and  hence  time  to  peak  in  a  hydrograph. 
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iv)  The  incorporation  of  headloss  due  to  expansion  and  contraction  of 
floodplain  flow  as  it  crosses  the  main  channel  is  best  achieved 
through  the  friction  headloss  computation. 

3 .1  Implications  for  the  improvement  of  MILHY 

The  conclusions  from  the  sensitivity  analysis  of  the  Ervine  and  Eliis  scheme 
isolate  friction  as  being  the  single  most  important  factor  in  the  prediction 
of  discharge  in  two-stage  channels.  Friction  is  identified,  therefore,  as 
being  the  key  to  improving  the  conveyance  capabilities  of  MILHY.  The 
analysis  showed  that  handling  of  frictional  headlosses  can  successfully 
incorporate  both  boundary  roughness  effects  and  the  effects  of  transverse 
currents  in  the  meandering  channels.  The  second  area,  and  worthy  of 
investigation  in  order  to  upgrade  the  predictive  performance  of  MILHY,  was 
the  impact  of  the  relatively  longer,  sinuous  path  length  of  the  main  channel 
over  the  floodplain. 

Three  key  processes  that  need  further  investigation  have,  therefore,  been 
identified;  these  are: 

1.  Improvement  of  the  handling  of  friction  to  incorporate  boundary 
roughness  and  transverse  circulations  within  the  main  channel 

2.  incorporation  of  turbulent  shear  stresses  between  the  main  channel 
and  floodplain  flow  segments 

3.  different  path  lengths  for  main  channel  and  floodplain  areas  to 
incorporate  sinuosity 

3.4.1  Incorporation  of  the  effects  of  turbulence  into  MILHY 

An  objective  identified  by  the  sensitivity  analysis  is  therefore  to 
incorporate  the  effects  of  turbulence  into  the  modelling  of  conveyance  in 
MILHY.  Two  significant  sources  of  turbulence  have  been  identified  as: 
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1)  apparent  shear  stresses  between  floodplain  and  channel  flow  segments 

2)  transverse  circulation  stresses  generated  at  meanders  o:  the  main 
channel 

The  relative  importance  of  these  two  sources  and  the  interaction  between 
them  is  not  clear  from  the  literature  although  it  is  dependent  on  the  depth 
of  flow  on  the  floodplain.  Because  the  analysis  of  the  Ervine  and  Ellis 
scheme  suggests  the  effects  of  transverse  circulations  at  meander  bends 
could  be  successfully  incorporated  in  modelling  of  boundary  friction,  it  was 
decided  to  concentrate  our  investigation  of  the  handling  of  turbulence 
modelling  on  the  shear  stresses  between  flow  segments.  This  investigation 
is  reported  in  Section  4.  It  is  hoped  that  the  Incorporation  of  the 
transverse  circulation  effects  can  be  Investigated  in  the  next  eighteen 
months. 

3.4.2  Incorporation  of  multiple  routing  pathways 

As  the  analysis  of  the  Ervine  and  Ellis  scheme  and  Fread  (1976)  suggest,  the 
different  path  lengths  of  the  sinuous  main  channel  and  straight  floodplain 
flows  will  affect  the  timing  of  a  floodwave  travelling  downstream.  At 
present,  MILHY  models  a  single  pathway  for  floodplain  and  channel  flows 
using  a  mean  channel  length  and  travel  timetable.  A  priority,  therefore, 
was  to  investigate  alternative  methods  of  incorporating  these  multiple 
pathways  of  flow  through  the  reach.  This  investigation  is  reported  in 
Section  5. 

3.5  Logical  development  of  a  research  scheme 

In  order  to  achieve  the  objectives  of  incorporating  the  two  processes 
identified  in  Section  3.4,  it  was  necessary  to  select  and  Implement 
processes  generally  only  modelled  in  hydraulically-based  schemes.  This 
raised  a  key  issue  in  the  MILHY3  project,  which  was  how  far  can  we 
successfully  model  channel  hydraulics  in  a  hydrological  model  and  when  does 
the  user  need  to  switch  to  a  hydraulically-based  model  such  as  RMA-2?  It 
seems  that  the  user's  decision  on  whether  to  use  a  hydrologically  or 
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hydraulically-based  scheme  has  been,  in  the  past,  based  paly  3n  his 
experience,  rather  than  on  a  clearly  defined  set  of  rules.  In  most 
applications,  the  user's  experience  is  sufficient  to  correctly  select 
hydrological  models  for  initial  passes  at  problems  and  ungauged 
applications,  and  hydraulic  models  for  detailed  engineering  applications. 
However,  there  remains  a  'middle-ground'  of  applications  where  data 
availability  may  be  difficult  or  its  quality  poor,  or  where  the  user 
requires  only  reach  outflow  data  but  for  a  complex  application. 

For  these  'middle-ground'  applications  it  is  not  clear  if  the  user  should 
pursue  a  hydraulically-based  scheme  with  often  extensive  set-up  and  run 
times,  or  would  be  better  selecting  a  simpler  hydrologically-based  scheme. 
The  margin  of  error  in  either  approach  for  these  'middle-ground' 
applications  needs  to  be  defined  so  that  a  series  of  operational  rules  may 
be  developed. 

The  logic  of  the  work  reported  here,  therefore,  is:- 

i)  having  identified  the  areas  worthy  of  investigation  in  order  to 
Improve  the  downstream  routing  capability  in  Section  3,  to 

11)  Implement  these  improvements  by  selecting  and  modifying  hydraulical 
processes  into  a  hydrologic  scheme  -  MILHY3 ,  then  to 

Hi)  test  and  validate  MILHY3  against  MILHY  and  MILHY2  on  the  River  Fulda 
catchment,  and  finally  to 

iv)  develop  a  set  of  operational  rules  for  the  application  of 

hydrological  and  hydraulically-based  schemes  by  a  comparative  study 
between  MILHY3  and  RMA-2. 

3.5.1  Implementation  of  research  scheme 

Having  achieved  the  point  aim  in  the  research  scheme  of  identifying  and 
selecting  areas  worthy  of  investigation,  the  next  stage  of  the  research  is 
to  Investigate  how  the  selected  two  processes,  momentum  transfer  between 
flow  segments  and  short-circuiting  of  floodplain  flows,  may  best  be 
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incorporated  into  MILHY.  These  two  processes  ire  investigated  individual lv 
in  Sections  4  and  5  where  alternative  approaches  are  examined,  the  tost 
appropriate  selected  and  a  series  of  initial  tests  run  to  check  the  validitv 
of  the  improvements  individually. 

Once  satisfied  that  the  incorporation  of  these  two  processes  improves  the 
predictive  capability  of  MILHY's  downstream  routing  component  when  applied 
individually,  the  performance  of  MILHY3  is  tested  as  a  rainfall-routing 
model  against  MILHY2  and  field  data  on  the  River  Fulda  catchment.  This 
analysis  is  reported  in  Section  6. 

Lastly,  the  comparative  study  between  MILHY3  and  RMA-2  is  reported  in 
Section  7. 


3.5.2  Work  identified  for  the  next  eighteen  months 

Section  3.4  indicated  that  there  are  other  important  processes  and 
parameters  worthy  of  further  investigation  that  are  beyond  the  scope  of  the 
three-year  research  period  reported  here.  The  key  processes  parameters 
identified  were  ranked  according  to  their  order  of  importance  and  the  two 
most  significant  processes  were  Implemented  and  validated.  The  investigation 
of  the  hydrology/hydraulic  model  trade-offs  was  not  anticipated  when  the 
original  research  proposal  was  submitted,  however  It  was  felt  that  this 
issue  was  of  Importance  for  the  MILHY  project  and  of  interest  in  may  other 
areas . 

It  is  hoped  in  the  next  eighteen  months  to  continue  the  hydrology/hydrauilc 
question  by  investigating  the  utility  in  linking  these  two  types  of  models. 
This  point  is  examined  further  in  Section  7.  In  the  next  eighteen  months  it 
will  also  be  possible  to  investigate  the  processes  active  in  an  in-bank 
meandering  channel. 
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4.  INCORPORATION'  OF  MOMENTUM  TRANSF-IR  3£TV;'£:.~ 

FLOODPLAIN-  AND  CHANNEL  SEGMENTS 

The  transfer  of  momentum  between  the  main  channel  and  floodplain  flow 
segments  was  identified  in  section  3.4  as  a  process  which  it  was  felt  would, 
if  incorporated  into  MILHY’s  downstream  routing  scheme,  make  a  significant 
improvement  to  the  overall  predictive  capability  of  MILHY.  The  objective  of 
the  work  reported  in  this  section,  therefore,  was  to  investigate,  implement 
and  validate  a  method  of  incorporating  momentum  transfer  between  flow 
segments  whilst  maintaining  MILHY's  parsimonious  data  requirements. 

4.1  The  hydraulics  of  momentum  transfer 


In  two-stage  channels,  the  irregular  cross-sectional  geometry  of  the  deep 
main  channel,  and  its  associated  shallow  floodplains,  generate  higher 
channel  velocities  in  the  main  channel  than  those  in  the  floodplain  flow 
segments.  This  is  due  to  the  relatively  greater  depth  of  flow  and  smaller 
wetted  perimeter  of  the  main  channel  in  comparison  to  the  floodplain. 

Figure  4.1  illustrates  the  velocity  isovels  for  a  two-stage  flame  experiment 
conducted  by  Knight  et  al.  (1983).  The  velocity  isovels  are  dimensionless 
parameters  because  the  observed  values  are  divided  by  the  mean  velocity  for 
cross-section,  where  V  »  Q/A.  The  difference  in  the  flow  velocity  between 
the  main  channel  and  floodplain  cause  a  transfer  of  longitudinal  momentum 
generally  from  the  main  channel  to  the  floodplain. 

The  physical  mechanisms  by  which  linear  momentum  is  transported 
perpendicular  to  the  direction  of  flew  are:- 

1)  secondary  currents 

2)  eddies  generated  In  the  mixing  zones  of  stream  tubes  of 

differing  veiocities 

3)  eddies  generated  by  flow  along  a  boundary 

4)  molecular  motion 

Wright  and  Carstens  (1970)  ranked  these  processes  on  a  scale  of  one  to  four 
In  order  of  their  effectiveness  at  transporting  momentum.  Figure  4.2  shows 
the  first  two  processes  listed  above,  and  Illustrates  the  position  of  the 
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Figure  4,2  :  Linear  momentum  transfer  by  secondary  currents  and 
turbulent  eddies  of  the  main  channel /f loodpla i n 
interface 
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eddies  generated  in  Che  mixing  zone  between  the  floodplain  and  the  main 
channel . 

As  suggested  earlier  in  section  3.1,  the  dominant  process  in  two-stage 
channels  is  the  turbulenc  eddy  in  the  mixing  zones  between  the  main  channel 
and  floodplain  (Pasche  and  Rouve,  1985).  Myers  (1978)  confirmed  this  when 
he  observed  chat  the  shear  stresses  obtained  in  two-stage  experimental  flume 
channels  were  much  greater  than  those  exerted  in  the  main  channel. 

The  momentum  transfer  distorts  the  shear  stress  profiles  on  the  beds  of  the 
main  channel  and  floodplain,  and  ir.ctc.cSes  the  shear  stresses  on  the 
floodplains  near  the  junction  with  the  main  channel.  Associated  with  the 
flow  interaction  are  high  shear  stresses  on  the  interfaces  between  the  main 
channel  and  the  floodplain  regions,  the  average  shear  stress  on  these 
interfaces  being  known  as  the  apparent  shear  stresses. 

A  great  deal  of  research  on  the  transfer  of  momentum  in  two-stage  channels 
has  been  carried  out  in  the  last  twenty-five  years,  starting  with  Sellin 
(1964)  and  Zheleznyakov  (1965).  Sellin  (1964)  was  the  first  to  identify  the 
turbulence  at  the  interface  between  main  channel  and  floodplain  by 
photographing  the  vortices  generated  by  the  turbulence  In  a  flume  study. 
Zheleznyakov  found  in  both  flume  (1965)  and  field  experiments  (1971)  that 
the  momentum  transfer  mechanism  decreased  the  overall  rate  of  discharge  for 
floodplain  depths  just  over  bankful.  Wright  and  Carstens  (1970)  described 
the  apparent  shear  stresses  as  a  drag,  retarding  main  channel  flows  and 
acting  as  a  propulsion  on  the  floodplain  flows.  Radojkovic  (1976)  pointed 
out  the  dependence  of  the  shear  stress  on  the  velocity  difference  between 
the  main  channel  whilst  Rajaratnam  and  Ahmadi  (1981)  concluded  that  If  the 
shear  velocity  were  known,  the  logarithmic  velocity  distribution  law  can  be 
applied  to  predict  floodplain  velocities  and  hence  discharges.  Recent 
research  into  the  behaviour  of  two-stage  channels,  for  example  Holden  and 
James  (1989),  has  concentrated  on  quantifying  the  physical  processes 
determining  the  rate  of  momentum  transfer  and  collecting  data  which  will 
allow  the  shear  stress  distribution  to  be  modelled.  In  the  United  Kingdom, 
the  Science  and  Engineering  Research  Council  are  funding  a  large  flume-based 
research  project  taking  place  In  four  universities  (Knight,  University  of 


/  t 


51 


52 


=  apparent  shear  stress  acting  upon  the  assumed  interface  i 
=  length  of  assumed  interface  i 
=  main  channel  solid  boundary  shear  force 

Figure  4.3  illustrates  these  forces  for  a  theoretical  example  where  the 
apparent  shear  stresses  are  assumed  to  be  acting  on  a  vertical  planar 
boundary  where  the  channel  and  floodplain  meet.  Rewriting  equations  4.1  and 
4.2  for  the  channel  segment  only  and  combining  them  with  equation  4.3  gives: 


where 
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In  any  application  A^,  P^  and  Sq  are  known  from  the  geometry  of  the 
cross-section  and  a  length  for  p^  can  be  assumed.  In  the  flume  the  average 
boundary  shear  stresses, may  be  measured  and  can  be  estimated 

using  equation  4.4.  The  discharge  for  each  flow  segment  can  then  be 
computed  from  the  corrected  retractive  forces. 

However,  in  an  ungauged  catchment,  it  is  extremely  unlikely  we  would  have 
boundary  shear  stress  data  available  and  unlikely  that  we  know  or  can 
estimate  the  length  of  the  apparent  shear  stress  boundary.  We  are  not 
suggesting,  therefore,  that  this  type  of  analysis  be  incorporated  into  MILHY 
but  investigation  of  the  application  of  this  method  in  flume  experiments 
does  provide  a  useful  Insight  into  the  relationship  between  the 
cross-sectional  geometry,  apparent  shear  stress  interfaces  and  accuracy  of 
the  discharge  prediction. 

4.2.2  Flume  experiments  investigating  apparent  shear  stresses 

Flume-based  research  programmes  provide,  at  present,  the  only  means  of 
collecting  data  on  the  distribution  boundary  shear  stresses  which  will 
enable  us  to  understand  and  later  model  the  processes  active  in  two-stage 
channels.  Field  data  of  two-stage  flood  events  are  notoriously  difficult 
and  sometimes  dangerous  to  collect.  The  variable  nature  of  flood  events 
means  that  flows  are  never  steady  enough  to  allow  even  reasonable 
measurement  of  the  velocity  fields,  while  field  measurement  of  boundary 
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shear  stresses  are  almost  impossible  to  collect.  Reliance  in  flume-based 
investigations  has  therefore  led  to  an  extensive  programme  of  modelling  a 
variety  of  geometrical  and  roughness  environments. 

In  these  flume  experiments  the  principal  objective  of  the  investigators  is 
to  develop  a  relationship  between  the  stage  and  discharge  in  the  main 
channel  and  floodplain  flow  segments.  The  investigators  hope  to  achieve 
this  objective  by  solving  equation  4.4  using  observed  flume  data  to  compute 
the  apparent  shear  stresses  from  the  solid  boundary  shear  stresses.  The 
boundary  shear  stresses  are  computed  using  either  the  Prandtl-von  Karman 
velocity  law,  utilizing  observed  velocity  data,  or  using  Patel's  (1965) 
relationship  between  head  difference  and  boundary  shear  stresses.  As  noted 
earlier,  however,  and  seen  in  equation  4.4,  the  value  of  the  computed 
apparent  shear  stress  is  dependent  on  the  length  of  the  assumed  interface 
over  which  the  apparent  shear  theoretically  acts.  Figure  4.3  illustrates  an 
assumed  Interface  in  a  vertical  plane,  a  method  which  has  been  utilized  by 
Chow  (1959)  and  Wright  and  Carstens  (1970).  Figure  4.4  shows  the  vertical 
plains  and  diagonal  interfaces  used  by  Womleaton  et  al.  (  1980),  and  Yen  and 
Overton  (1973).  and  the  horizontal  Interfaces  utilized  by  Deuller  et  al. 
(1967). 

Wormleaton  et  al.  (1982)  carried  out  a  comparative  investigation  of  the 
apparent  shear  stresses  computed  over  the  three  types  of  planar  interface. 
Their  results  were  reported  as  an  apparent  shear  stress  ratio,  that  is  the 
ratio  of  the  apparent  shear  stress  to  the  average  shear  stress  including  the 
assumed  interface.  As  the  apparent  shear  stress  tends  to  zero,  the  ratio 
will  also  tend  to  zero,  implying  no  shear  on  the  interface.  The  results  of 
Wormleaton  et  al.  (1982)  are  shown  on  Figure  4.5  a,  b  and  c  for  the 
vertical,  diagonal  and  horizontal  interfaces  respectively,  where  the 
apparent  shear  stress  ratio  and  an  Inundation  ratio  are  compared.  The 
inundation  ratio  is  defined  as  the  depth  of  flow  on  the  floodplain  divided 
by  the  depth  of  the  main  channel.  The  series  A,  B,  C  and  D  illustrate  the 
effects  of  Increasing  the  floodplain  Manning's  roughness  from  0.011  for 
series  A  through  0.014  (B),  0.017  (C),  to  0.021  series  D. 

Analysis  of  Figure  4.5  a,  b  and  c  shows  that  the  apparent  shear  stress 
declines  with  increasing  depth  of  flow  on  the  floodplain  in  all  three  planar 
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ngure  4.4  :  /ertioal,  diagonal  ana  horizontal  apparent 
shear  stress  interfaces 


Figure  4.6 


Angle  inclination  of  zero  shear  stress 
interfaces 
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interfaces.  The  order  of  magnit.de  difference  between  the  apparent  shear 
stresses  computed  for  the  vertical  interfaces  and  those  computed  on  the 
diagonal  and  horizontal  interfaces  should  also  be  noted.  This  shows  that 
the  vertical  interface  is  much  nearer  to  the  turbulent  eddies  photographed 
by  Sellin  (1964).  Analysis  of  the  boundary  stress  distributions  showed  that 
the  negative  apparent  shear  stress  ratios  computed  for  the  diagonal  and 
horizontal  interfaces  at  higher  floodplain  inundation  depths  indicate  a 
transfer  of  momentum  from  the  zone  >f  flow  above  the  main  channel  to  the 
within-bank  main  channel  zone. 

Wormleaton  et  al.  (1982)  wary  of  the  criticism  that  all  their  apparent  shear 
stress  values  were  computed  using  a  single  cross-sectional  width,  developed 
a  relationship  by  regression  analysis  between  geometric  and  velocity 
parameters  for  the  apparent  shear  stress.  This  could  then  be  compared  with 
data  collected  by  other  authors  often  for  very  different  applications  and  so 
utilize  data  from  a  wide  variety  of  cross-sectional  geometries.  Wormleaton 
et  al.  (1982)  give  a  final  regression  equation  for  the  form  of  a  vertical 


interface  as: 


where  A  V  is  the  velocity  difference  between  the  floodplain  and  main  channel 
flow  segments,  computed  from  the  Manning  equation  (equation  3.19). 

Utilizing  the  data  from  34  experimental  set-ups  the  coefficient  of 
determination  for  equation  4.5  was  0.983.  Data  collected  by  Myers  (1978), 
Crory  and  Elsawy  (1980)  and  Ghosh  and  Jena  (1971)  were  found  to  conform 
closely  with  relationship  4.5. 

Yen  and  Overton  (1973)  tackled  the  problem  from  an  alternative  perspective, 
using  the  measured  boundary  9hear  stress  profiles  to  position  an  interface 
along  which  no  shear  would  take  place.  The  cross-section  could  then  be 
divided  up  using  these  no-shear  boundaries  and  the  discharge  computed  easily 
as  it  would  be  directly  related  to  the  segment's  cross-sectional  area.  Yen 
and  Overton  (1973)  attempted  to  relate  the  angle  of  a  zero  shear  Interface, 
pivoting  around  the  main  channel/floodplain  intercept  (see  Figure  4.5)  to 
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observed  discharge  values.  If  this  angle  could  then  be  related  to  the 
cross-sectional  geometric  parameters,  then  this  method  could  be  applied  verv 
simply  to  a  wide  variety  of  problems. 

Yen  and  Overton's  (1973)  results  showed  that  the  angle  of  inclination  o:  the 
zero  shear  stress  plane  varied  with  both  the  ratio  of  floodplain  to  main 
channel  width,  and  the  ratio  of  floodplain  inundation  to  main  channel  depth. 
With  a  range  of  width  ratios  between  2.2  and  5.4  the  angle  of  inclination 
varied  by  as  much  as  20°,  with  the  angle  increasing  as  the  width  ratio 
decreased.  The  angle  of  inclination  varied  with  a  depth  ratio  range  from 
0.2  to  1.8  by  60°  with  the  angle  increasing  linearly  from  15°  to  50°,  with 
the  depth  ratio  up  to  a  depth  ratio  of  approximately  1.0  when  the 
relationship  becomes  exponential.  The  angle  of  inclination  of  zero  shear 
stress  for  a  particular  cross-section  does  not  vary,  therefore,  when  the 
depth  ratio  is  above  2. 

The  results  of  Wormleaton  et  al  (1982)  reported  in  Figure  4.6  agree  with 
those  of  Yen  and  Overton  (1973)  and  show,  therefore,  that  when  the  ratio  of 
the  floodplain  inundation  to  main  channel  depth  is  approximately  2  or  above, 
the  two-stage  channel  may  be  considered  as  a  single  system.  Below  this 
ratio  the  distribution  of  the  turbulent  shear  stresses  has  been  shown  to  be 
complex  where  no  one  single  position  of  the  apparent  shear  stress  interface 
or  stress  ratios  can  be  adequately  applied  to  describe  the  boundary  shear 
stresses  over  a  variety  of  cross-sectional  geometries. 

4.2.3  Implications  of  flume-based  experiments  for  the  prediction 
of  the  discharge  capacity  of  two-stage  channels 

It  was  noted  earlier  that  the  main  reason  that  the  relationship  between 
cross-sectional  area  and  discharge  does  not  hold  for  two-stage  is  the 
transfer  of  momentum  between  the  main  channel  and  the  floodplain.  The 
flume-based  experiments  reported  in  Section  4.2.2  attempt  to  quantify  these 
momentum  transfers  by  balancing  the  gravitational  and  retarding  forces  by 
the  Introduction  of  an  apparent  shear  stress  over  a  dividing  interface 
between  segments  of  flow.  However,  in  order  to  compute  the  discharge 
capacity,  we  need  to  develop  a  relationship  between  easily  measured 
geometric  parameters  and  the  stage/discharge  rating  curve.  There  are 
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several  alternatives  that  can  be  used:- 

1)  We  could  use  empirical  relationships  developed  from  flume  experiments 
to  predict  the  percentages  of  flow  in  each  cross-sectional  segment. 
These  are  developed  from  regression  analysis  of  the  computed  apparent 
shear  stresses  on  assumed  interfaces.  Examples  include  the 
relationships  developed  by  Wormleaton  et  al.  (1982)  (reported 
earlier,  see  equation  4.5)  and  by  Knight  and  Demetriou  (1983). 

2)  We  could  attempt  to  divide  the  cross-section  using  the  zero-shear 
interfaces,  suggested  by  Yen  and  Overton  (1973). 

3)  We  could  divide  the  cross-section  using  shear  interfaces  and  make 
some  assumption  about  the  amount  of  momentum  transfer  across  these 
interfaces . 

Each  of  these  alternatives  are  now  considered.  The  first  proposition  to  use 
empirically  developed  relationships  seems  attractive  in  that  it  would  be 
simple  to  apply.  However,  the  relationships  have  been  developed  using  data 
collected  in  flume  experiments  which  have  had  limited  cross-sectional 
geometries.  Table  4.1  shows  the  geometric  parameters  of  the  major  flume 
investigations  that  have  published  this  type  of  data.  Comparison  of  the 
floodplain  Co  main  channel  widths  shows  a  maximum  ratio  of  3  where  in  the 
River  Fulda  catchment,  flood  Inundation  maps  Illustrate  a  ratio  of  up  to  50. 
Similarly  the  maximum  Manning's  roughness  applied  to  the  floodplain  Is 
0.022,  whilst  Chow  (1959)  suggests  a  typical  grazed  pasture  to  have  a 
Manning's  'n'  value  of  0.03.  To  generate  empirical  relationships  applicable 
to  the  sorts  of  two-stage  channels  typical  in  Europe,  therefore,  there  is  a 
need  for  further  flume  experiments  with  much  wider  and  rougher  floodplains. 
Until  this  Is  achieved,  it  would  be  Inadvisable  to  extrapolate  the  existing 
relationships  to  geometries  and  roughness  outside  those  reported  in  Table 
4.1. 


The  second  alternative  given  is  to  divide  the  cross-section  along  zero-shear 
Interfaces,  as  suggested  by  Yen  and  Overton  (1973).  As  there  is  no  momentum 
transfer  across  the  zero-shear  Interface,  the  Manning  equation  will  hold  for 


/_ 
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each  cross-sectional  flow  segment.  Although  Yen  and  Overton  computed  the 
angle  of  Incidence  of  the  interfaces  (see  Figure  1.6)  for  width  ratios  up  to 
5,  the  sensicivity  of  this  angle  to  floodplain  roughness,  means  the  results 
cannot  be  reliably  applied.  Computing  the  area  of  flow  segments  based  on 
the  angle  of  inclination  around  the  main  channel/ floodplain  interface  is 
also  rather  more  difficult  than  if  a  vertical,  horizontal  or  diagonal 
interface  could  be  used. 

However,  it  is  important  to  note  here  that  the  zero  shear  interface  has  been 
applied  widely  for  a  number  of  years,  as  Lotter’s  (1933)  technique,  where 
the  cross-section  is  divided  vertically  into  a  number  of  segments,  and  where 
the  interfaces  are  ignored  in  computing  the  wetted  perimeter  of  an 
interface.  As  Yen  and  Overton  (1973)  have  shown,  though,  such  interfaces 
are  not  vertical,  but  move  from  the  inclined  towards  the  horizontal  as  the 
depth  of  flow  increases.  Zero  shear  interfaces  can  be  applied,  therefore, 
for  vertical,  diagonal  and  horizontal  inclinations  by  ignoring  the  assumed 
interface  in  the  wetted  perimeter  computation  and  taking  the  solid 
boundaries  only. 

The  third  suggestion  to  compute  the  discharge  capacity  of  the  cross-section 
was  to  divide  the  cross-section  using  the  shear  interfaces,  making  an 
assumption  about  the  amount  of  momentum  transfer  across  these  interfaces. 

In  a  similar  way  to  the  zero  shear  interfaces,  shear  interfaces  have  been 
applied  in  a  great  number  of  environments,  using  vertical,  diagonal  and 
horizontal  inclinations.  The  assumption  here  is  that  the  apparent  shear 
stress  Is  equal  to  the  average  shear  stress  (i.e.  that  of  Wormleaton  et  al. 
(1982),  apparent  shear  stress  ratio  is  equal  to  1,  see  Figure  4.5),  so  that 
Che  interface  can  be  included  as  part  of  the  wetted  perimeter  in  the 
discharge  capacity  computation. 

Wormeaton  et  al.  (1982)  computed  the  discharge  for  the  zero-shear  and  shear 
interfaces  for  ail  three  inclinations  over  a  variety  of  floodplain 
roughnesses  up  to  n  ■  0.021.  Their  results  showed  that,  as  expected,  the 
computed  discharge  values  converged  to,  or  were  smaller  than,  the  observed 
values  when  the  f loodplain/channel  depth  ratio  increased  to  2,  for  all 
Interface  inclinations.  However,  the  accuracy  of  the  discharge  prediction 
using  these  six  techniques  was  considered  only  with  variation  in  the  depth 
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ratio  and  floodplain  roughness;  the  width  ratios  were  not  considered. 

The  implication  of  the  flume-based  experiments  to  the  computation  of 
discharge  in  two-stage  channels,  is  that  no  one  single  technique  of 
incorporating  turbulent  exchange  between  the  main  channel  and  floodplain,  is 
appropriate  for  all  geometric  and  roughness  environments.  The 
flume  experiments  need  to  be  extended  to  represent  field  environments  with 
wider  and  rougher  floodplains  before  a  set  of  operational  rules  on  the 
suitability  of  zero-shear  or  shear  interfaces  and  their  angle  of 
inclination,  can  be  developed. 

4 .3  Incorporation  of  momentum  transfer  into  MILHY 


Analysis  of  the  flume-based  experiments,  in  section  4,2,  has  shown  that 
there  is  no  single  method  of  incorporating  momentum  transfer  between  flow 
segments  that  is  appropriate  for  all  cross-sectional  geometries  and 
roughnesses.  For  this  reason,  and  because  of  the  lack  of  comparative  work 
on  wide  and  rough  floodplains,  it  was  decided  to  incorporate  a  number  of 
different  methods  into  MILHY  and  test  the  accuracy  of  the  discharge 
predictions  against  observed  field  data  collected  from  the  River  Fulda 
catchment . 

4.3.1  Selection  of  methods  for  incorporation  into  MILHY 

Four  methods  of  dividing  the  cross-section  to  incorporate  momentum  transfer 
were  considered.  These  were:- 

1.  Vertical  subdivision,  with  zero  shear  interfaces. 

2.  Vertical  subdivision,  with  an  apparent  shear  stress  ratio  »  1. 

3.  Diagonal  subdivision,  with  zero  shear  interfaces. 

4.  Diagonal  subdivision,  with  an  apparent  shear  stress  ratio  *  1. 

At  present,  method  2,  that  is  vertlvcal  subdivision  with  an  apparent  shear 
stress  ratio  equal  to  l,  is  Incorporated  into  both  MILHY  and  MILHY2 .  By 
application  of  these  four  different  methods  it  should  be  possible  to  test 
the  sensitivity  of  the  generated  rating-curve  to  the  interface  Inclination 
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and  apparent  shear  stress  ratio.  This  sensitivity  could  then  be  compared  n 
the  impact  on  the  rating  curve  of  variation  in  the  cross-sectional  geomecrv 
and  roughness  parameters.  If  the  analysis  showed  the  rating-curve  to  be 
sensitive  to  the  computational  method,  then  further  methods  including 
horizontally  inclined  and  Yen  and  Overton's  (1973)  angle  of  Inclination 
could  be  incorporated  and  tested. 

4.3.2  Incorporation  of  the  four  methods  into  MILHY 

The  four  methods,  identified  above,  of  incorporating  momentum  transfer 
between  the  main  channel  and  floodplain,  are  the  same  four  methods  utilized 
by  Knight  and  Hamed  (1984),  Knight  and  Hamed  (1984)  tested  the  accuracy  of 
Che  four  identified  techniques  in  predicting  discharge  by  comparing  the 
predicted  results  with  those  collected  in  flume  experiments  conducted  by 
Knight  and  Demetriou  (1983),  reported  in  Table  4.1.  For  consistency,  then, 
and  to  ensure  the  correct  cross-sectional  definitions  were  being  applied  to 
MILHY  for  each  of  the  four  methods,  the  equations  of  definitions  reported  in 
Knight  and  Hamed' s  (1984)  paper  were  incorporated  into  MILHY.  These 
equations  are  given  in  Table  4.2  whilst  Figure  4.7  defines  the 
cross-sectional  geometry  variables  used.  Analysis  of  the  equations  in  Table 

4.2  shows  that  the  wetted  perimeter  of  the  interface  is  included  In  the  main 
channel  computation  in  methods  2  and  4,  whilst  being  excluded  in  methods  1 
and  3 . 

These  four  methods  were  then  incorporated  Into  the  rating  curve  generation 
routine  (subroutine  CMPRC) ,  Introducing  an  option  variable  into  the  'datal' 
dataset.  These  changes  are  recorded  in  the  source  code  of  MILHY3  which  is 
given  in  Appendix  3.  It  is  noted  that  the  cross-section  definitions 
reported  in  Table  4.2  are  only  incorporated  for  stage  elevations  above 
bankful.  This  is  significant,  especially  for  method  4,  where  the  wetted 
perimeter  of  the  main  channel  would  be  otherwise  effectively  extended  above 
the  stage  level. 

4.4  Sensitivity  of  the  rating  curve  to  Interface  inclination 


There  were  several  objectives  In  undertaking  a  sensitivity  analysis  of  the 
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Table  4.2 


Alternative  geometric  definitions  to  incorporate  segment  interactions 
(after  Knight  and  Hamed ,  1984) 


Method 

Flood  Plain 

Main  Channel 

Area 

Wetted  Perimeter 

Area  Wetted  Perimeter 

1 

(H-h) 

(B-b) 

B-b  +  H-h 

2bH 

2b  +  2h 

2 

(H-h) 

(B-b) 

B-b  +  2( H-h) 

2bH 

2b  +  2H 

3 

(H-h) 

( B-b/2) 

B-b  +  H-h 

b(  B+h) 

2b  +  2h 

4 

(H-h) 

(B-b/2) 

B-b  +  H-h 

b( H+h) 

2b  +  2h  + 

->  o  i  / 

2((H-h)‘>b“)  '  2 

racing  curve.  these  were: 


1)  co  escabLish  whether  any  one  method  improved  Che  accuracy  of 

the  racing  curve  in  comparison  to  observed  field  rating  curves, 
for  a  field  cross-section 


2)  to  establish  whether  there  is  a  significant  difference  in  the 
predicted  rating  curves  generated  by  each  of  the  four  methods  for 
wide  floodplains  with  greater  boundary  roughnesses  than  those 
reported  in  Table  4.1. 

3)  to  compare  the  difference  in  the  computed  rating  curve  attributable 
to  the  interface  inclination  method,  with  the  difference  due  to 
variability  in  the  cross-sectional  geometry  and  roughness  parameters 

To  answer  these  three  questions,  it  was  necessary  to  apply  the  four 
interface  inclination  methods  to  both  field  cross-sections,  to  achieve 
objective  one,  and  hypothetical  reaches,  to  achieve  objectives  two  and 
three.  Whilst  the  field  cross-sections  are  similar  to  the  theoretical 
cross-sections  as  they  have  broadly  rectangular  main  channels  and  flat  wide 
floodplains  (see  Figure  4.9)  application  of  field  cross-sections  provided 
the  only  comparison  to  an  observed  rating  curve  possible.  Objectives  two 
and  three  can  be  achieved  by  comparison  of  the  divergence  in  predicted 
rating-curve  discharge  between  computation  methods  applied  to  hypothetical 
cross-sections . 


plication  of  the  four  Interface  inclination  methods 


The  cross-section  at  Bad  Hersfeld  on  the  River  Fulda,  West  Germany,  was 
selected  in  order  to  compare  the  accuracy  of  the  four  computation  methods 
against  a  field  rating  curve.  The  rating  curve  at  Bad  Hersfeld  was  extended 
to  out-or-bank  conditions  using  data  from  gauged  extreme  events  for 
floodplain  inundation  depths  of  up  to  3.2  metres.  This  depth  corresponds 
approximately  co  the  1  In  100  year  event.  At  Bad  Hersfeld  the  floodplains 
are  symmetrical  about  the  main  channel  with  a  floodplain  to  main  channel 
width  ratio  (B/b)  of  10.  The  bankful  depth  (h)  is  4.1  m  whilst  the 
floodplains  on  either  side  of  the  main  channel  are  pasture  grazed  with 


Figure  4.9  :  Cross-sect,ion  on  River  Fulda  at  Meeklar,  half  way 
between  Bad  Hersfeld  and  Rotenburg 
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cattle.  The  four  interface  inclination  methods  we  re  applied  fir  three  seta 
of  geometric  and  roughness  environments  and  the  discharge  at  increments  of 
3.5  m  computed.  The  rating  curve  was  also  computed  for  the  first  two  cases 
with  the  cross-section  being  treated  as  a  single  system,  chat  is  with  no 
interfaces  to  divide  the  cross-section  into  segments.  The  rating  curves 
produced  from  these  three  applications  are  reported  in  Tables  .  3  to  -.5. 

A  theoretical  cross-section  was  established  to  achieve  objectives  two  and 
three  noted  above,  with  a  rectangular  main  channel  and  a  floodplain  rise 
from  channel  to  valley  side  of  only  0.1  metres.  The  floodplain  to  main 
width  ratios  considered  were  10  and  20,  as  noted  earlier.  Flame  experiments 
by  numerous  authors  have  investigated  smaller  width  ratios.  Wormleaton  et 
al.  (1982)  reported  that  discharge  predictions  from  all  the  interface 
inclination  methods,  vertical,  diagonal  and  horizontal,  converged  on  a 
common  solution  as  the  floodplain  inundation  depth  to  main  channel  depth 
(H/h)  approached  2.  To  check,  this,  discharge  predictions  were  calculated 
for  depth  ratios  (H/h)  up  to  2.2  were  computed  at  0.5m  stage  increments.  As 
well  as  the  four  interface  inclination  methods,  the  rating  curve  was 
computed  treating  the  cross-section  as  a  single  flow  segment.  Where 
friction  or  slope  parameters  varied  between  main  channel  and  floodplain 
segments,  a  mean  average  between  the  two  values  was  applied  to  the  single 
segment  case.  This  was  true  for  both  the  hypothetical  and  Bad  Hersfeld 
cross-sections.  The  hypothetical  cross-section  results  are  reported  in 
Tables  4.6  to  4.10. 

4.4.2  Sensitivity  of  the  rating  curve  to  the  computational  stage  increment 

Whilst  the  sensitivity  analysis  for  the  interface  inclination  methods  was 
being  set  up,  it  was  noticed  that  the  discharge  predictions  were  sensitive 
to  the  computational  stage  increment.  The  rating  curve  is  computed  at 
twenty  evenly  spaced  elevation  Increments,  these  increments  being  generated 
from  the  minimum  and  maximum  cross-sectional  elevations.  Figure  4.8 
Illustrates  the  Impact  of  relatively  small  changes  In  the  computational 
stage  Increment,  from  0.29  to  0.61m,  on  the  discharge  values  predicted  for 
the  Bad  Hersfeld  cross-section.  The  maximum  difference  in  the  predicted 
discharges  between  stage  increment  values  occur  at  1.5m  above  bankful  where 
the  discrepancy  is  approximately  150  m^s  *.  The  observed  rating  curve  at 
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this  elevation  gives  a  discharge  value  of  250  m^s-".  Figure  1.3  also  shows 
Chat  a  decrease  in  the  stage  increment  size  does  not  necessarily  improve  the 
accuracy,  as  the  rating  curves  for  increment  sizes  of  2.29m  and  ".him  are 
almost  identical. 

The  maximum  elevation  in  the  field  cross-sections  of  the  River  Fulda 
catchment  are  determined  by  the  valley  side.  The  computational  stage 
increment  therefore  is  computed  from  the  height  of  the  constraining  side 
(valley)  wall. 

4.4.2  Results  of  the  sensitivity  analysis 

The  results  of  the  sensitivity  analysis  of  the  computed  rating  curve  to  the 
interface  inclination  and  variation  in  geometric  parameters,  are  tabulated 
In  Tables  4.3  to  4.10.  Tables  4.3  to  4,5  show  the  results  computed  for  the 
Bad  Hersfeld  cross-section  and  also  record  the  percentage  error  of  each  of 
the  interface  Inclination  methods  against  an  observed  rating  curve.  Tables 
4.6  to  4.10  show  the  results  for  a  hypothetical  cross-section,  and  the 
percentage  error  in  these  tables  Indicate  the  deviation  from  the  MILHY 
solution  as  no  'observed'  rating  curve  was  available. 

Table  4.3  contains  the  observed  discharge  values  and  the  computed  values 
from  MILHY2  and  interface  inclination  methods  1  to  4.  Table  4.3  confirms 
that  MILKY2  Incorporates  method  2,  and  In  further  tables,  therefore,  both 
are  not  shown.  Manning's  'n'  values  of  0.035  for  the  main  channel  and 
floodplain  were  selected  for  the  first  run  reported  in  Table  4.3.  This 
value  corresponds  to  the  tabulated  values  suggested  in  Chow  (1959).  The 
channel  and  floodplain  slopes  were  set  at  0.0006,  computed  from  the  field 
rating  curves  from  Bad  Hersfeld  and  Rotenburg,  the  next  gauging  station 
downstream. 

Results  from  the  Bad  Hersfeld  station 

Table  4.3  shows  the  discharge  predictions  from  the  four  interface  methods 
computed  using  the  parameter  values  reported  above.  The  mean  average  error 
of  the  discharge  predictions  over  the  observed  figures  was  computed  for  each 
method  over  a  range  of  inundation  depths. 
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Table  4.3  shows  that  all  nee  hods  i  -  ill  inundation  depths  overpredicced  the 
carrying  capacity  of  the  cross-section.  The  average  error  shows  that  tie t hod 
2,  the  method  utilized  by  MILHY2,  gave  the  worst  prediction.  The  best 
overall  prediction  was  given  by  the  single  segment  method,  however  this  was 
not  so  surprising  as  both  the  boundary  roughness  and  slope  variables  were 
constant  across  Che  section. 

Table  4.3  also  shows  that  there  was  no  consistent  difference  In  predictive 
performance  between  the  methods  incorporating  the  shear  face,  methods  2  and 
4,  and  the  zero  shear  methods.  Also  worthy  of  note  Is  that  the  percentage 
error  increases  with  depth  in  all  methods  except  method  3.  It  Is  important 
Co  remember  here  that  the  algorithm  Co  reduce  Manning's  n  as  depth  increases 
has  been  removed  from  MILHY3,  for  reasons  specified  in  section  3.4.1. 

Comparison  of  Tables  4.3  and  4.4  shows  how  increasing  the  floodplain 
boundary  roughness  can  more  than  half  the  error  of  the  predictions  for  all 
methods.  The  difference  in  mean  average  errors  between  computation  methods 
is,  however,  the  same  as  those  in  Table  4.3.  This  suggests  that  the 
carrying  capacity  computation  is  more  sensitive  to  the  boundary  roughness 
value  than  the  form  of  Che  main  channel/ f loodplain  interface. 

Table  4.4  also  shows  that  the  percentage  error  does  not  increase  with 
Increasing  floodplain  inundation  depth,  as  suggested  by  Table  4.3.  In  Table 
4.4  the  percentage  errors  values  indicate  that  all  four  computation  methods 
are  converging  to  the  observed  discharge  as  the  inundation  depth  increases 
and  approaches  the  main  channel  depth.  This  suggests  that  a  floodplain 
boundary  roughness  value  of  0.07  corresponds  more  closely  to  the  field 
conditions  than  the  initial  value  used  by  0.035.  The  logic  behind  this 
argument  lies  in  chat  as  floodplain  inundation  depth  Increases  to  the  main 
channel  depth,  the  two-stage  channel  behaves  as  a  single  system  and 
therefore  all  of  the  computation  methods  should  converge  on  a  common 
solution.  If  they  do  not,  as  in  Table  4.3,  this  suggests  that  the  initial 
variable  values  used  are  not  realistic. 

Table  4.5  shows  the  effects  of  incorporating  meandering  in  the  channel  by 
reducing  the  slope  value  used  to  0.0001  from  0.0006.  This  value  is 
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calculated  from  the  ratio  of  Che  main  channel  length  to  Che  /allev  length 
between  Bad  Hersfeld  and  Rotenburg  on  the  River  Fulda.  Comparison  of  Tables 

4.3  and  4.5  shows  that  reducing  the  slope  of  the  main  channel  improves  the 
predictions  of  the  carrying  capacity  in  all  computation  methods.  However, 
we  are  not  looking  to  calibrate  this  particular  applications,  rather  we  are 
aiming  to  compare  the  efects  of  variation  in  the  parameters  against  the 
method  of  computation.  The  results  of  this  field  application  suggest  that 
computation  method  utilized  in  MILHY2  (i.e.  method  2)  give  in  Tables  -.3  and 

4.4  the  poorest  prediction  of  the  carrying  capacity  of  the  cross-sect  ion . 
However,  Tables  4.4  and  4.5  show  that  the  prediction  can  be  improved  to  much 
greater  extent  by  more  accurate  selection  of  parameter  values  than  by 
altering  the  computation  technique. 

Results  from  a  hypothetical  cross-section 

Tables  4.6  to  4.10  report  the  predictions  of  the  carrying  capacity  for  a 
hypothetical  cross-section,  comparing  methods  one  to  four  and  discharges 
computed  by  treating  the  cross-section  as  a  single  segment.  The  percentage 
error  values  reported  are  computed  from  the  MILHY2  predictions,  which 
utilizes  method  2.  The  percentage  error  values  allow  comparison  of  the 
relative  sensitivity  of  the  discharge  predictions  to  variation  in  the 
computation  method  and  parameters.  The  absolute  accuracy  of  the  techniques 
cannot  be  computed  as  this  is  a  hypothetical  application.  It  is  not  useful, 
therefore,  to  directly  compare  the  percentage  errors  from  the  Bad  Hersfeld 
section  to  the  hypothetical  application. 

Analysis  of  Tables  4.6  to  4.10  shows  that  method  1  produces  a  very  close 
approximation  to  the  predictions  produced  from  the  MILHY2  computations , 
under  all  of  the  boundary  roughness  and  geometry  environments.  In  all 
cases,  methods  3  and  4  rank  second  and  third  respectively  in  their  closeness 
to  the  MILHY2  predictions.  Methods  1,  3  and  4  under-predict  the  carrying 
capacity  In  comparison  to  the  MILHY2  predictions  in  all  five  examples. 
Comparison  of  Tables  4.6  and  4.10  where  the  f loodplaln/main  channel  width 
ratio  has  been  increased  to  20  from  10,  shows  that  this  increase  in  the 
width  ratio  has  made  little  impact  on  the  comparative  accuracy  of  the 
computation  methods.  There  has  been  no  radical  change  in  the  difference  in 
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the  mean  average  errors  between  the  four  computation  methods. 

Comparison  of  the  hypothetical  and  3ad  Hersfeld  applications 

Analysis  of  the  two  sets  of  results  has  shown  that  the  method  incorporated 
into  MILHY2,  method  2,  gives  the  largest  prediction  of  the  carrying  capacitv 
of  the  cross-section  in  both  the  Bad  Hersfeld  and  hypothetical  sections. 

The  3ad  Hersfeld  section  results  suggest  that  method  2  gives  the  worst 
prediction  of  the  four  methods,  which  all  over-predict  the  carrying 
capacity.  This  suggests  that  all  four  methods  do  not  introduce  enough 
friction  over  the  assumed  interfaces  between  the  main  channel  and  floodplain 
to  mimic  the  retarding  effects  of  momentum  exchange.  Method  4  assumes  a 
diagonal  interface  and  an  apparent  shear  stress  ratio  equal  to  one,  and 
introduces  the  most  additional  boundary  friction  of  the  methods,  hence 
producing  the  lowest  prediction  of  carrying  capacity  (see  Tables  4.3  to 
4.10).  This  suggests  that  in  the  field  apparent  shear  stress  ratios  on 
diagonal  interfaces  may  be  greater  than  one,  rather  than  less  than  one  as 
Jonnleaton  et  al.  (  1982)  found  in  Fig.  4.5b.  Alternatively,  these  results 
suggest  that  the  true  position  of  the  interface  is  between  the  vertical  and 
diagonal,  as  apparent  shear  stress  ratios  on  the  vertical  interface  are  very 
much  greater  than  one  (see  Figure  4.5a).  Apparent  shear  stress  ratios  of 
greater  than  one  could  be  incorporated  into  the  MILHY  scheme  by  multiplying 
the  wetted  perimeter  of  the  apparent  interface  in  the  main  hannel 
computation  until  the  ratio  was  reduced  to  one. 

4.4.4  Conclusions 


From  the  analysis  of  the  results  above,  it  is  possible  to  make  several 
conclusions : 

1)  The  three  methods  utilized  to  incorporate  turbulent  exchange  between 
the  main  channel  and  floodplain,  more  accurately  predict  the  carrying 
capacity  of  a  cross-section  than  the  technique  used  in  MILHY2. 

2)  .All  four  methods  over-predicted  the  carrying  capacity  because  they 
failed  to  Introduce  enough  additional  boundary  friction  to  mimic  the 
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effaces  of  turbulent  exchange .  Method  -  Introduce!  the  tost 
additional  friction  and  therefore  gave  the  best  predictions. 

3)  Increasing  the  boundary  roughness,  Manning's  'n*  for  the  flood  pi  sins 

was  more  effective  at  reducing  the  over-prediction  of  the  carrying 
capacity  of  the  section  than  increasing  the  wetted  perimeter  of  the 
interface,  or  assuming  an  apparent  shear  stress  ratio  of  one. 

i . 5  Implications  for  the  improvement  of  MILHY 

The  initial  sensitivity  analysis,  reported  above,  of  the  alternative  methods 
of  incorporating  turbulent  exchange,  has  shown  that  these  simple 
redefinitions  can  improve  the  predictive  capability  of  MILHY.  Application 
of  the  techniques  to  wider  and  rougher  floodplains  than  previously  reported 
(see  Table  4.1)  and  to  cross-sectional  field  geometries  which  3re  only 
broadly  rectangular  has  been  successful.  Therefore,  these  turbulent 
exchange  routines  are  worthy  of  inclusion  for  validation  of  Che  whole  MILHY3 
scheme,  ,4s,  however,  these  routines  have  only  been  tested  against  one  field 
section  in  this  analysis,  it  was  felt  that  all  the  routines  should  be 
carried  onto  the  final  analysis,  and  not  just  the  best  method,  method  4, 
identified  in  this  analysis.  Although  it  is  not  our  intencion  to  prove  with 
any  statistical  certainty  the  most  appropria  e  technique,  as  this  would 
equLre  possible  hundreds  of  applications,  we  do  intend  to  be  able  to 
suggest  guidelines  for  these  techniques.  We  anticipate  that  a  combination 
of  turbulent  exchange  routines  and  a  multiple  routing  routine  (reported  in 
Section  5)  may  highlight  a  different  exchange  routine  than  the  one  selected 
here.  It  is  also  possible  that  in  a  rainfall,  runoff  and  route  application 
Che  turbulent  exchange  routine  invoked  may  become  unimportant.  For  these 
reasons,  then,  all  four  of  the  turbulent  exchange  routines  are  included  in 
the  sensitivity  analysis  reported  in  Section  6. 


33 


SECTION'  5  :  I NCORPORAT 1 3 S'  IF  ••g'LTI?-.;  RAFTING  REACHES 


The  sensitivity  analysis  of  Ervine  and  Ellis  (1937),  reported  in  Section 
3,  found  that  Che  sinuosicy  of  Che  min  channel  was  important  in 
determining  Che  length  of  the  flow  path.  In  two-stage  channels, 
therefore,  multiple  routing  reaches  provide  a  means  of  incorpora t ir.g  the 
differing  path  lengths  of  the  sinuous  main  channel  and  relatively 
straight  floodplain. 

5 . 1  The  behaviour  of  two-stage  channel  flow 

In  two-stage  channels  there  is  a  tendency  for  floodplain  flow  to 
"short-circuit"  the  general  mote  sinuous  route  taken  by  the  main  channel 
(Fread,  1976).  The  infrequency  of  out-of-bank  flows  means  that  flows 
taking  the  shortest  path  downslope  do  not  develop  the  secondary  flow 
system  necessary  for  the  development  of  meanders.  Einstein  and  Shen 
(1964)  suggested  that  the  secondary  flow  system  is  initiated  itself  by 
shear  along  a  rough  bank.  The  shorter  path  lengtn  of  the  floodplain  flow 
is  exacerbated  by  the  steeper  gradient  of  the  floodplain  in  comparison 
with  the  main  channel,  giving  faster  velocities  and  travel  rimes  for 
floodwaves  passing  downstream. 

The  accelerating  effects  of  the  path  length  and  slope  on  floodplain 
fLows  are  diminished,  however,  by  effects  of  boundary  friction.  If 
floodplain  flow  depths  are  small  then  the  hydraulic  radius  will  also  be 
small  and  hence  velocities  reduced.  Floodplain  boundary  roughnesses  also 
tend  to  be  higher  than  those  in  the  main  channel  because  of  vegetation 
and  obstructions  such  as  hedges. 

5.1.1  Comparison  of  main  channel  and  floodplain  boundary  roughness 

As  noted  earlier  in  Section  3.1.3  the  retarding  effects  of  boundary 
roughness  tend  to  decline  as  the  hydraulic  radius  or  stage  Increases. 
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This  is  particularly  true  for  the  broadly  rectangular  -tain  : ha n t. e  1  s  found 
in  the  River  Fulda  catchment.  On  the  floodplains,  however,  the  situation 
is  complicated  by  vegetation  and  nan-nade  structures.  As  the  floodplain 
inundation  depth  increases  debris  can  become  trapped  in  hedges  and  fences 
and  the  boundary  roughnesses  may  increase.  Klaassen  and  dwaari  (19'-.) 
showed  that  the  spacing  of  hedges  and  trees  is  critical  in  computing  the 
friction  of  floodplains. 

When  selecting  Mannings  ’n'  values  for  the  floodplain  segments  in  the 
rating  curve  computation  of  MILHY,  it  is  essential  to  consider  not  only 
the  general  land  use  but  also  the  spacing  and  height  of  any  hedges  or 
fences . 


5.1.2  Objective  of  investigating  multiple  routing  reaches 

The  aim  of  incorporating  multiple  routing  reaches  is  to  separate  the 
conflicting  effects  of  the  straighter  pathway  and  higher  boundary 
friction  of  floodplain  flow,  on  the  selection  of  the  most  appropriate 
Mannings  'n'  value.  By  removing  considerations  of  the  sinuosity  of 
floodplain  flows,  the  selection  of  the  correct  'n'  value  should  be 
simplified.  As  the  sensitivity  of  MILHY3  to  the  Manning  'n'  value  is 
anticipated,  by  improving  the  selection  of  a  suitable  value  and  improving 
the  representation  of  the  physical  processes  active  in  the  two-stage 
channel,  the  predictive  capability  of  MILHY  sboulld  be  improved. 

5.2  Modelling  alternatives 

The  optimal  method  of  modelling  the  conveyance  of  a  floodplain  through  a 
two-stage  channel  would  be  to  use  a  two  or  even  possibly 
three-dimensional  finite  element  model.  Such  models  allow  for  turbulent 
exchange  between  elemental  areas  and  so  predict  a  pattern  of  flow  across 
the  section.  The  amalgamation  of  these  complex  and  time-consuming  finite 
element  models  with  simple  hydrologic  models,  such  as  MILHY,  is  discussed 
further  In  Section  7. 
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The  objective  in  this  section  is  to  develop  a  simple  jne-J  imer.s i  ona  1 
technique  that  can  be  incorporated  into  MILHY.  Such  a  technique  should 
therefore  have  - 

little  additional  data  requirements 
low  computer  processor  demands 
be  capable  of  validation 


These  demands  leave  several  alternative  approaches  available.  These 
are 


1)  To  develop  a  stage/reach  length  relationship.  This  approach  was 

suggested  by  Perkins  (1970),  when  he  incorporated  into  a 
ratnfall/runof f  model  a  routine  to  increase  reach  length  linearly 
from  the  main  channel  thalweg  distance  at  bankful  to  the  shortest 
reach  length  dictated  by  the  floodplain  slope,  at  the  maximum 
stage. 


2)  To  develop  an  empirical  adjustment  to  the  roughness  coefficients 

of  the  floodpLain  and  main  channel.  This  approach  was  suggested 
by  Tlngsanchali  and  Ackermann  (1976),  where  Manning's  'n'  value 
was  weighted  by  a  ratio  of  reach  lengths  between  the  actual 
floodplain  distance  and  a  schematized  straight  floodplain  and  main 
channel.  Such  that 


n* 


5.1 


where  n*  =  adjusted  Manning's  'n' 

n^  =  Manning's  'n'  floodplain 

Lj  ■  reach  length  of  floodplain 

L  =  reach  length  of  main  channel 
me  ° 


I 


3) 


Replace  the  variable  storage  coefficient  routing  technique  in 
MILHY,  with  a  St.  Venant  technique  utilizing  a  weighted  four-point 
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implicit  finite  difference  solution,  modified  by  Freud  (  ;  376 )  ti 
incorporate  the  differing  path  lengths  of  fi.odpLain  r.J  main 
channel  flows. 

-i)  Separate  floodplain  and  main  channel  flows  and  route  jsing  the 

existing  routines  in  MILHY  the  flow  downstream  assuming  no 
exchange  of  flow  along  the  reach. 

The  simplest  solution  to  apply  is  approach  four,  where  each  cross- 
sectional  segment,  used  to  develop  the  rating  curve,  is  touted 
individually  downstream.  There  are  several  disadvantages  however 

1)  flow  has  to  be  apportioned  to  floodplain  or  channel  as  the 

top  of  the  reach,  and  these  proportions  are  fixed  throughout 
the  reach.  This  implies  the  assumption  that  the  cross-sectional 
geometry  is  fairly  constant  downstream; 

it)  there  is  no  exchange  of  momentum  between  the  main  channel  and 
floodplain  along  the  reach; 

ill)  floodplain  flows  cannot  cross  main  channel  flows. 

Despite  these  disadvantages,  approach  four  seemed  to  be  a  logical  first 
step  Into  tackling  the  problem  of  floodplain  flows  "short-circuiting"  the 
main  channel.  Exchange  of  momentum  between  the  main  channel  and 
floodplain  has  been  Incorporated  at  the  valley-sections  (see  Section  i) , 
and  it  was  felt  Important  at  this  stage  to  compare  the  sensitivity  of  the 
outflow  hydrograph  to  the  effects  of  variability  In  the  turbulent 
exchange  routines  or  the  multiple  routing  of  floodplain  and  channel 
flows.  If  the  downstream  "short-circuiting"  effects  were  identified  as 
being  significant,  then  it  would  be  appropriate  to  investigate  Perkins 
approach  as  another  simple  alternative,  or  a  more  radical  replacement  of 
the  routing  subroutine  with  Fread  (1976)  St.  Venant  solution. 
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5.3 


Incorporation  of  multiple  routing  into  MILKY 


The  incorporation  of  multiple  routing  into  MILKY  proved  to  be  a 
relative!/  simple  matter,  with  changes  in  the  source  code  being  required 
mainly  to  facilitate  the  apportioning  of  flow  across  the  cross-section. 
Each  time  step  of  the  inflow  hydrograph  was  apportioned  into  the 
cross-sectional  segment  according  to  the  stage  it  reached  at  the 
cross-section  nearest  the  top  of  the  reach.  It  was  necessary,  therefore, 
to  develop  stage-rating  curves  and  compute  the  percentage  of  total  flow 
with  stage  for  each  segment.  These  computations  are  included  in  the 
CMPRC  subroutine  and  are  output  to  the  results  file.  The  inflow 
hydrograph  is  then  apportioned  in  the  ROUTE  subroutine  when  the 
individual  segment  rating  curves  and  percentage  curves  are  recalled.  For 
each  segment,  a  TRAVEL  TIME  and  ROUTE  command  are  invoiced  and  the  outflow 
hydrographs  are  then  added  to  give  the  total  discharge  across  the 
segment . 

Tliis  approach  minimised  the  changes  required  in  the  source  code  as  most 
of  the  multiple  routing  can  be  achieved  by  repetition  of  commands  in  the 
'datal'  dataset. 

5.4  Application  of  multiple  routing  reaches 

It  was  important  at  this  stage  to  test  the  impact  of  multiple  routing  on 
the  outflow  hydrograph.  The  relative  importance  of  the  technique 
compared  to  other  routine  modifications  and  the  sensitivity  of  the  whole 
scheme  to  parameter  variability  is  investigated  in  Section  6. 

Multiple  routing  was  applied,  therefore,  to  a  theoretical  reach  with 
rectangular  cross-sectional  geometry  assumed  to  be  constant  downstream, 
and  a  reach  from  the  River  Fulda,  between  Bad  Hersfeld  and  Rotenburg.  A 
variety  of  inflow  hydrographs  were  applied  to  the  theoretical  reach,  in 
order  to  look  at  the  impact  of  the  depth  of  inundation  on  the  travel  time 
of  the  floodplain  and  the  effects  on  Che  outflow  hydrograph.  The  River 
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Fulda,  however,  provided  fieLd  data  against  which  various  rjaghr.ess  me 
routing  lengths  could  be  tested  but  with  a  Limited  number  of  observed 
flood  in  flow  and  outflow  hydrographs.  A  L  in  19  year  event  was 
available  and  enough  flood  frequency  data  was  available  to  generate  a 
in  190  year  event,  assuming  a  similar  shape  as  the  1  in  1 '  year  event. 

The  1  in  190  year  event  could  Legitimately  be  used  to  compare  the 
accuracy  of  the  predicted  peak  stage  and  discharge. 

All  the  simulations  reported  involve  the  routing  of  an  input  hydrograph 
from  between  two  cross-sectional  stations.  The  runoEf  contribution  of 
the  drainage  area  between  the  two  stations  is  not  considered.  The 
results  of  these  simulations  are  reported  in  Tables  5.1  to  5.3  and 
Figures  5.2  to  5.7.  The  tables  show  the  time  to  peak,  peak  discharge  and 
maximum  floodplain  inundation  of  the  outflow  (downstream)  hydrograph. 

5.4.1  Application  to  the  3ad  Hersfeld  to  Rotenburg  reach 

The  results  from  the  application  of  multiple  routing  reaches  to  the  Bad 
Hersfeld  to  Rotenburg  reach  are  found  in  Tables  5.1  and  5.2,  and  Figures 
5.2  to  5.5.  As  reported  earlier  in  Section  4.4.1,  the  cross-sectional 
geometry  at  Bad  Hersfeld  is  broadly  rectangular  with  the  floodpLains 
being  symmetrical  about  the  main  channel.  The  reach  from  Bad  Hersfeld  to 
Rotenburg  is  approximately  24  km  (15  miles)  in  length  with  a  sinuous  main 
channel;  this  can  be  seen  on  Figure  5.1.  At  Rotenburg  the  bankfull  depth 
is  4.8  m  as  compared  tv  4.1  m  at  Bad  Hersfeld,  with  a  bankfull  discharge 
of  180  m^s  The  valley  section  Is  asymmetrical  at  Rotenburg  with  the 
left  hand  floodplain  being  approximately  300  metres  wide  whilst  the  right 
hand  floodplain  rises  steeply.  The  bankfull  width  at  Rotenburg  is 
approximately  50  m  as  compared  with  30  m  at  Bad  Hersfeld.  When  multiple 
routing  Is  Invoked,  therefore,  the  observed  hydrograph  at  3ad  Hersfeld  Is 
apportioned  to  floodplain  and  main  channel  segments  according  to  the 
rating  curve  developed  for  the  Bad  Hersfeld  cross-section.  The  travel 
time  table  is  then  developed  for  each  cross-sectional  segment  using 
whichever  of  the  two  segments  produces  the  smaller  rating  curve.  The 
maximum  floodplain  inundation  values  reported  in  Tables  5.1  and  5...  are 
computed  at  the  downstream  end  of  the  reach,  that  is  at  Rotenburg. 
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j  Table  5.1 

l 

;  Bad  Hersfeld  to  Rotenburg  reach,  1  In  10  year  event 


Time  to  Peak 
( hours) 

Peak  Discharge 
(m»  ) 

Maximum  Floodplain 
Inundation 

m 

Observed 

38 

407 

0.33 

MILHY2 

38 

285 

0.09 

Multiple  routing 

40 

330 

0.17 

Multiple  routing 
floodplain  length  5% 

40 

333 

0.18 

Multiple  routing 
floodplain  length  30% 

40 

352 

0.21 

Multiple  routing 
floodplain  'n'  30% 

38 

355 

0.22 
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Table  5.2 


I 


Bad  Hersfeld  to  Rotenbarg  reach,  1  in  100  year  event 


Time  to  Peak 
(hours) 

Peak  Discharge 

(mV1) 

Maximum  Floodplain 
Inundation 

m 

Observed 

38 

744 

0.90 

MILHY2 

38 

665 

0.78 

Multiple  routing 

40 

634 

0.73 

Multiple  routing 
floodplain  length 

36 

3  0Z 

684 

0.81 

Multiple  routing 
floodplain  'n'  305! 

36 

668 

0.78 
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the  1  la  10  year  event  is  shown  in  Figure  5.2,  indicating  a  travel  tine 
of  the  peak  discharge  of  approximately  nine  hours.  The  inflow  hydrograph 
at  Bad  Hersfeld  has  been  scaled  up,  in  line  with  the  flood  frequency  data 
available  to  provide  tne  1  in  100  year  event  and  consequently  the  1  in 
100  year  event  has  the  same  form  as  the  1  in  10  year  event.  At  Bad 
Hersfeld  the  1  in  100  year  event  corresponds  to  an  increase  in  the 
floodplain  inundation  depth  of  approximately  1  m  over  the  1  in  10  year 
event . 

Figure  5.3  compares  the  observed  outflow  hydrograph  at  Rotenburg  with  the 
outflow  hydrographs  produced  by  MILHY2  and  the  multiple  routing 
technique.  The  greatest  difference  in  the  three  hydrographs  occurs  in  the 
over-bank  section  which  is  the  area  of  particular  interest.  The 
corresponding  time  to  peak,  peak  discharge  and  maximum  inundation  depths 
of  these  three  hydrographs  are  recorded  on  Table  5.1.  Both  Che  figure 
and  table  show  that  the  single  routing  technique  used  in  MILHY2 
effectively  smooths  the  inflow  hydrograph  to  too  great  an  extent.  This 
reduces  the  peak  discharge  and  floodplain  inundation  depth.  The  effects 
of  multiple  routing  are  to  sharpen  up  the  hydrograph,  thereby  Increasing 
Che  peak  discharge  and  inundation  depth.  The  multiple  routing  technique 
halves  the  MILHY2  errors  in  both  the  peak  discharge  and  inundation  depth. 
Table  5.1  also  shows  that  the  multiple  routing  technique  produced  a  time 
to  peak  of  40  hours,  two  hours  later  than  the  observed  peak.  However, 
the  observed  Inflow  and  outflow  hydrographs  were  digitised  at  three 
hourly  Intervals  and  therefore  errors  of  less  than  three  hours  <_an  be 
effectively  Ignored. 

As  reported  earlier  in  this  section,  the  main  objective  of  incorporating 
multiple  routing  reaches  was  to  simulate  the  effects  of  the 
short-circuiting  of  floodplain  flow,  reducing  the  floodplain  reach 
length.  In  the  next  simulation  reported  in  Table  5.1,  therefore,  the 
reach  length  of  the  floodplain  segments  was  reduced  by  5%.  This  produced 
only  very  small  changes  in  the  outflow  hydrograph  in  comparison  with  Che 
multiple  routing  hydrograph  shown  on  Figure  5.4.  Analysis  of  the  flood 
Inundation  maps  available  for  the  River  Fulda  indicated,  however,  that 
the  floodplain  reach  length  may  be  up  to  30%  shorter  than  the  main 
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channel.  The  hydrograph  produced  by  reducing  the  floodplain  reach  length 
by  30%  is  shown  on  Figure  5.4,  indicated  by  the  triangles.  Comparison  of 
Figures  5.3  and  5.4  shows  that  reducing  the  floodplain  length  by  30% 
makes  a  significant  improvement  on  the  prediction.  Figure  5.4  shows  as 
well,  however,  that  a  similar  effect  can  be  achieved  by  reducing  the 
Mannings  'n'  roughness  coefficient  by  30%.  Chow  (1959)  showed  that  the 
effects  of  sinuosity  of  a  channel  can  alter  the  'n'  value  by  up  to  30%. 

As  noted  earlier,  however,  the  aim  of  this  investigation  is  to 
incorporate  the  processes  active  in  two-stage  channels  and  remove  the 
reliance  on  empirical  coefficients.  The  effectiveness  of  these  new 
process  'modules'  in  improving  the  predictive  capabilities  of  MILHY  in 
comparison  with  the  established  techniques  used  is  investigated  in 
Section  6. 

Table  5.2  and  Figure  5.5  show  the  simulation  results  for  the  1  in  100 
year  event  on  the  River  Fulda  reach.  In  contrast  to  the  1  in  10  year 
event  the  MILHY2  prediction  gives  a  higher  peak  discharge  result  than  the 
multiple  routing  reach.  The  proportional  difference  between  the  MILHY2 
and  multiple  routing  technique  is,  however,  much  smaller  in  the  1  in  100 
year  storm  being  4%,  whereas  the  1  in  10  year  difference  was  11%.  This 
suggests  that  as  the  floodplain  inundation  depth  increases  the 
cross-section  behaves  as  a  single  system.  However,  predicting  the 
floodplain  reach  length  or  Mannings  'n'  value  has  the  same  effects  on  the 
1  in  100  year  event  as  the  1  in  10  year  event,  the  peak  discharge  is 
increased . 

5.4.2  Application  to  a  hypothetical  reach 


The  aim  of  investigating  the  impact  of  multiple  routing  on  a  hypothetical 
reach  was  to  look  at  the  relative  impact  of  the  floodplain  inundation 
depth  on  the  outflow  hydrograph.  A  hypothetical  reach  was  set  up  with 
symmetrical  rectangular  cross-sections  at  upstream  and  downstream 
stations  with  a  floodplain/main  channel  width  ratio  of  10.  The  main 
channel  was  a  constant  depth  of  2.4  metres,  so  that  the  main  channel 
capacity  remained  constant  downstream.  This  meant  that  the  proportion  of 
flow  on  the  floodplain  was  correct  throughout  the  reach  and,  therefore, 
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the  analysis  could  concentrate  on  the  effects  of  inundation  depth. 

Table  5.3  and  Figures  5.6  and  5.7  summarize  the  results  from  this 
investigation  into  the  impact  of  floodplain  inundation  on  depth  on  the 
outflow  hydrograph.  In  all  these  simulations  the  floodplain  and  channel 
reach  length  were  held  constant  at  20  km.  The  seven  inflow  hydrographs 
were  generated  by  scaling  the  1  in  10  year  observed  hydrograph  from  Bad 
Hersfeld  (see  Figure  5.2).  The  scaling  factor  used  for  each  storm  is 
recorded  in  column  two  of  Table  5.3.  By  using  the  same  form  of  inflow 
hydrograph,  it  was  hoped  to  be  able  to  isolate  the  effects  of  the 
inundation  depth  on  the  time  to  peak  of  the  outflow  hydrograph. 

The  results  in  Table  5.3  are  summarized  on  Figure  5.8  which  plots  the 
percentage  error  between  the  MILHY2  and  multiple  routing  predictions  of 
peak  discharge.  Negative  errors  occur  when  the  multiple  routing 
predictions  are  less  than  the  MILHY2  prediction,  positive  errors  occur 
when  the  multiple  routing  predictions  are  greater.  Figure  5.9  plots  the 
error  between  the  MILHY2  and  multiple  routing  technique  when  the 
floodplain  routing  reach  length  is  reduced  by  30%.  Analysis  of  Figures 
5.8  and  5.9  show  that  the  predictions  from  the  two  techniques  converge  as 
the  floodplain/ main  channel  depth  ratio  increases  to  1  (equal  to  a 
Worraleaton  et  al.  (198z)  depth  ratio  of  2).  The  maximum  error  between 
the  two  techniques  occurs  when  the  floodplain  inundation  depth  is  0.3  of 
the  bankfull  depth. 

5.4.3  Conclusions 


1.  The  maximum  impact  of  the  multiple  routing  technique  occurs 
when  floodplain  inundation  depths  are  small. 

2.  At  these  small  Inundation  depths  (depth  ratios  0.3)  the 
utilization  of  multiple  routing  significantly  improves 
the  prediction  of  the  peak  discharge,  errors  are  halved. 

3.  Reducing  the  floodplain  routing  length  by  30%  reduced 
the  travel  time  of  the  peak  uiscuarge. 


a. 
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Table  5.3 

Hypothetical  reach 


Storm 

Multiple 

Time  to  peak, 
(hours) 

Peak  discharge 

(mV1) 

Maximum 

floodplain 

inundation 

m 

1 

1 

MILHY2 

36 

309.8 

1.16 

M.  Routing 

36 

317.4 

1.18 

9 

L 

0.1 

MILHY2 

36 

34.8 

- 

M.  Routing 

36 

32.8 

- 

3 

0.2 

MILHY2 

42 

57.2 

0.09 

M.  Routing 

36 

55.0 

0.07 

4 

0.5 

MILHY2 

40 

139.5 

0.61 

M.  Routing 

38 

126.9 

0.56 

5 

1.5 

MILHY2 

36 

536.0 

1.68 

M.  Routing 

36 

540.8 

1.74 

6 

2 

MILHY2 

34 

682.0 

1.98 

M.  Routing 

34 

685.0 

1.99 

7 

3 

MILHY2 

34 

1061.0 

2.65 

M.  Routing 


34 


1063.0 


2.65 


TIME  (hours) 


(hours) 


FLOODPLAIN/MAIN  CHANNEL  DEPTH  RATIO 


FLOODPLAIN/MAIN  CHANNEL  DEPTH  RATIO 
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SECTION  6  :  VALIDATION  OF  MILHY3 


Before  MILHY3  can  be  deemed  to  be  operational  the  following  questions 
must  be  answered: 

1.  Do  the  mathematical  algorithms  introduced  represent  the  processes 
we  are  trying  to  model? 

2.  Are  the  mathematical  algorithms  robust? 

3.  Is  the  accuracy  of  the  predicted  outflow  hydrograph  a  significant 
improvement  over  earlier  versions  of  MILHY? 

4.  Is  the  resolution  of  each  new  submodel  consistent  with  each  other, 
with  existing  submodels  and  appropriate  for  ungauged 
applications? 

5.  Can  a  set  of  operational  rules  be  developed  for  MILHY3? 

The  behaviour  of  the  results  reported  in  the  initial  analysis  in  Sections 
4  and  5  suggests  that  the  turbulent  exchange  and  multiple  routing 
algorithms  are  correctly  modelling  the  effects  of  friction  on  the  outflow 
hydrograph.  These  simulations  prove  that  the  algorithms  are  internally 
valid,  that  is  that  the  outflow  hydrograph  does  not  change  when  the  input 
parameters  are  held  constant  (Hermann,  1967),  and  that  the  algorithms  are 
mathematically  robust.  Simulations  also  showed,  although  not  reported 
here,  that  the  coding  of  the  turbulent  exchange  and  multiple  routing 
routines  was  correct,  in  that  during  in-bank  conditions  neither  routine 
was  invoked.  Sections  4  and  5  have  therefore  satisfactorily  verified  the 
coding  of  the  new  submodels  in  answer  to  questions  1  and  2,  and  show  that 
the  new  routines  are  modelling  the  processes  they  were  designed  to  and 
are  robust.  The  results  also  suggest  that  the  new  submodels  are  an 
improvement  on  the  predictions  made  by  earlier  MILHY  versions,  although 
further  testing  is  required  to  thoroughly  establish  this. 
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In  the  selection  of  the  new  algorithms  the  objective  was  to  improve  the 
predictive  capability  of  MILHY  by  incorporating  the  representation  of  the 
active  processes.  However,  the  selection  was  made  taking  into 
consideration  the  limitations  of  the  final  objective  of  maintaining  MILHY 
as  an  operational,  ungauged  model.  These  limitations  were: 

(1)  the  data  requirements  of  the  new  algorithms  should  be  small; 
in  particular  field  work  should  not  be  required 

(2)  any  additional  demands  made  of  the  user  in  the  establishment  of 
the  datasets  should  not  require  detailed  hydrological  knowledge 
of  the  physical  processes  or  computer  expertise 

(3)  the  computer  demands,  in  terms  of  CPU  and  operating  space,  of  the 
new  routines  should  not  limit  the  models  portability  to  the 
IBM-PC  level 

In  answer  to  question  4,  posed  above,  these  limitations  reduce  the  risk 
of  over-development  in  the  resolution  of  modelling  in  a  particular 
process  area.  The  concentration  of  model  resolution  on  the  most 
important  processes  in  an  otherwise  lumped  model,  improves  efficiency  but 
has  inherent  dangers.  The  dangers  include  limiting  the  portability  of  a 
model  by  increasing  the  quality  and  quantity  of  its  data  demands.  It  is 
felt  therefore  that  the  submodels  developed  in  sections  4  and  5  are  of  a 
consistent  resolution  with  each  other,  the  existing  submodels  and  the 
limitations  of  ungauged  modelling. 

Before  a  set  of  operational  rules  (question  5)  can  be  developed,  a  series 
of  applications  must  be  made.  It  is  to  this  question  and  to  the 
quantification  of  the  possible  improvement  in  the  predictive  accuracy  of 
MILHY  that  the  rest  of  this  section  is  directed.  In  the  analysis  of  the 
applications  that  follow,  all  five  of  these  questions  will  be 
investigated  and  the  conclusions  drawn  from  the  results  of  sections  4  and 
5  reassessed. 

To  answer  the  two  outstanding  questions  posed  at  the  start  of  this 


section,  several  specific  questions  may  be  asked.  These  are:- 

(1)  Is  the  outflow  hydrograph  more  sensitive  to  variability  in  its 
parameters,  or  to  the  process  submodels  utilized? 

(2)  What  is  the  relative  impact  of  the  submodels  introduced  in  MILHY3 
in  comparison  with  the  infiltration  algorithm  introduced  in 
MILHY2? 

(3)  What  is  the  impact  on  the  outflow  hydrograph  of  the  conflicting 
effects  of  the  new  submodels? 

(4)  What  is  the  effect  of  the  scale  of  the  catchment  on  the  three 
questions  posed  above? 

To  answer  these  and  to  complete  the  answering  of  the  general  questions 
asked  earlier,  it  was  felt  necessary  to  undertake  a  sensitivity  analysis 
of  MILHY3.  This  analysis  would  compare  hydrograph  predictions  made  by 
MILHY3  and  MILHY2  against  data  collected  from  a  gauged  catchment  for 
out-of-bank  conditions. 

It  is  Important  to  note  at  this  point  that  the  infiltration  algorithm 
introduced  in  MILHY2  has  only  been  applied  to  single  subcatchments.  Work 
by  Anderson  and  Howes  (1986)  reported  in  Section  1  of  this  report,  showed 
that  the  infiltration  algorithm  significantly  improved  the  generation  of 
the  runoff  hydrograph.  The  relative  importance  of  this  improvement  in  a 
large  catchment,  where  several  subcatchraents  are  utilized,  has  not  yet 
been  tested. 

6.1  Selection  of  a  field  site 


From  the  objectives  listed  above,  a  list  of  prerequisites  for  the  study 
catchment  was  developed. 

6.1.1  Prerequisites  of  a  study  catchment 


1. 


The  catchment  must  be  in  a  temperate  region  with  a  minimum  of 
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forested  areas. 

The  catchment  must  not  exceed  a  maximum  area  of  approximately 
2500  km"  (1000  sq.mi.),  (Williams,  1975). 

There  should  be  a  minimum  of  man-made  interferences  with  the 
hydrology  of  the  catchment,  such  as  reservoirs  or  land  drainage 
schemes . 

There  should  be  intermediate  gauging  stations  throughout  the 
catchment,  so  that  a  suite  of  catchments  can  be  developed. 

The  catchment  must  be  subject  to  floodplain  inundation. 

The  floodplains  should  be  geometrically  simple,  bounded  and 
preferably  with  a  similar  land-use  type  throughout  the  catchment. 

Data  should  be  easily  available  and  reliable.  The  minimum 
requirements  are:- 

i)  topographic  maps 

ii)  soils  classification  maps 

iii)  cross-sectional  geometries  for  a  number  of  gauging 
stations  in  the  catchment 

iv)  rating  curves  for  these  gauging  stations 

v)  a  number  of  storm  hydrographs  including  out-of-bank 
events  from  each  gauging  station 

vi)  corresponding  rainfall  data  for  these  storm  events  and 
the  week  preceding.  Data  must  consist  of  a  minimum  of 
daily  Information  from  at  least  one  rain  gauge  In  the 
catchment 

vtl)  flood  frequency  data  to  indicate  the  frequency  and  size  of 
out-of-bank  events 
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6.1.2  The  River  Fulda  Catchment 


These  seven  prerequisites  limited  prospective  study  catchments  to  the 
rural  regions  of  Western  Europe  and  areas  of  the  U.S.A.  From  a  short 
list  of  regions  meeting  the  prerequisites,  the  River  Fulda  catchment  in 
West  Germany  (see  Figure  6.1)  was  selected,  primarily  because  of  the 
efficiency  and  rapid  response  to  requests  made  to  the  relevant  water 
authorities  and  meteorological  offices.  As  well  as  the  prerequisite 
data,  the  local  authorities  in  the  River  Fulda  catchment  were  able  to 
provide  :- 


i)  an  outline  of  the  extent  of  floodplain  inundation  for  a  storm 
event  in  1946,  corresponding  to  the  1  in  200  year  event 

ii)  daily  precipitation  values  for  approximately  45  rain  gauge 
stations  (see  Figure  6.2) 

ill)  continuous  rainfalL  data  for  two  stations,  Sad  Hersfeld  and 
Kunzel 1-Die tershausen 

iv)  for  one  storm,  the  water-equivalent  of  snow,  daily  minimum  and 
maximum  temperatures,  relative  humidities  and  cloud  cover 

v)  long-profiles  of  two  of  the  reaches,  between  Bad  Hersfeld  and 
Rotenburg,  on  the  River  Fulda,  and  between  Marbach  and 
Hermanns p legal ,  on  the  River  Haune  (a  tributary  of  the  River 
Fulda) 


At  this  poLnt,  we  would  like  to  acknowledge  the  help  and  cooperation  of 
the  Water  Authority,  Wasserwirtschaf tsamt ,  Fulda  for  the  provision  of  the 
hydrological  data  and  the  Meteorological  Office,  Deutcher  Wetterdienst 
Eentralarat,  Offenbach,  Frankfurt,  for  the  meteorological  records, 
collected  during  three  visits  to  the  catchment  in  the  period  November 
1986  to  June  1988.  A  copy  of  this  raw  data  has  been  forwarded  to  the 
Environmental  Laboratory,  WES,  Vicksburg,  Mississippi,  who  kindly 
provided  the  soils  classification  maps. 
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6 .2  Establishment  of  the  River  Fulda  Catchment 

The  River  Fulda  catchment  to  Rotenburg  consists  of  a  drainage  area  of 
7 

approximately  2523  km  (974  sq.mi.)>  drained  by  the  River  Fulda  and  its 
tributaries.  The  main  tributary  is  the  River  Haune  which  joins  the  Fulda 
at  Bad  Hersfeld;  in  addition,  the  River  Luder  joins  the  Fulda  at  Lutterz. 
There  are  eight  river  gauging  stations  in  the  catchments,  marked  on 
Figure  6.1  for  which  six  storm  events  have  been  collated.  The  positions 
of  the  gauging  stations  have  enabled  the  division  of  the  catchment  into 
nine  subcatchments  depicted  in  Figure  6.3. 

During  the  visits  to  the  catchment,  sketches  were  made  and  photographs 
taken  that  enabled  the  technical  channel  cross-sections  to  be  extended 
across  the  floodplains.  Estimates  were  also  made  during  these  visits  of 
the  Mannings  'n'  roughness  values  of  the  channel  and  floodplain 
throughout  the  catchment.  Figures  6.4  to  6.6  are  photographs  taken  at 
Hetterhausen,  Unter-Schwarz  and  Rotenburg,  and  show  the  topography  and 
land-uses  typical  throughout  the  catchment.  The  photographs  show  that:- 

i)  in  the  upper  reaches  the  channel  is  tree-lined 

ii)  the  floodplains  are  extensive  and  flat 

iii)  the  floodplains  throughout  the  catchment  are  vegetated  by  short 
grass 

iv)  there  are  few  obstructions  on  the  floodplains,  there  are  few 

fences,  and  the  small  villages  tend  to  have  been  built  clear  of 
the  areas  subject  to  flooding 

v)  the  channel  is  broadly  rectangular  in  cross-section 

vi)  the  channel  is  sinuous 

Tables  6.1  and  6.2  collate  some  of  the  topographic  dimensions  of  the 
subcatchments  and  the  channel  geometries  at  the  gauging  stations. 
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Table  6.1 


Parameter  Values  for 

Subcatchments  in  the 

River  Fulda 

Catchment 

K 

\ 

Subcatchment  Area 

km2 

Max.  elev. 

m . 

Mia.  elev. 

m . 

Main  channel 

length 

l 

\ 

km . 

1 

401  56 

838 

365 

14 

402  506 

550 

232 

36 

40  3 

182 

700 

232 

25 

404 

469 

775 

216 

n  7 

405 

394 

416 

193 

33 

406 

148 

700 

265 

■>  u 

407 

274 

610 

209 

3-4 

408 

90 

518 

193 

9 

409 

403 

391 

179 

24 

Iota  1 

2523 

838 

179 

22  7 

I 

i 

I 


l _ L 
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Table  6.2 


Parameter  values  for  Gauging  Stations  in  the  River  Fulda  Catchment 


Station 

Bankful  1 

depth 

m. 

Bankf ull 

width 

m. 

Bankful 

capacity 

m3S-1 

Hetterhausen 

2.3 

17.0 

26 

Karamerzell 

2.0 

20 . 1 

33 

Lutterz 

3.2 

18.0 

18 

Unter-Schwarz 

3.0 

18.0 

50 

Marbach 

2.3 

8.0 

10 

Hermannsplegal 

2.5 

16.5 

22 

Bad  Hersfeld 

4.1 

30.3 

76 

Rotenburg 

4.8 

50.0 

179 

Note  : 


Cross-sectional  geometrical  data 
Unter-Schwarz  gauging  stations, 
estimates  made  during  visits  to 


was  not  available  for  the 
The  figures  that  appear  are 
the  catchment. 


The  six  storm  events  were  identified  as  being  discrete  events,  that  is 
where  the  hydrographs  rose  and  fell  back  to  baseflow  condition  with  a 
single  seven  day  record.  For  each  of  these  events,  the  daily  rainfall 
totals  for  the  three  preceding  weeks  was  collected  in  order  to  compute 
the  antecedent  conditions. 

In  order  to  compute  the  rainfall  in  each  of  the  nine  subcatchments,  the 
theisson  polygon  technique  was  used  to  weight  the  daily  rainfall  total 
from  each  of  the  45  rain  gauges  shown  in  Figure  6.2.  Polygons  of  the 
area  that  could  be  associated  with  a  particular  raingauge  were  drawn  as 
if  the  catchment  were  a  peneplain. 

Table  6.3  shows  the  percentage  occurrence  of  each  of  the  major  soil 
groups  in  each  of  the  nine  sub-catchments.  A  certain  amount  of 
interpolation  and  generalization  occurred  during  the  computation  of  this 
table,  as  the  pixel  definition  of  the  soils  classification  map  (1  pixel  = 
100  metres)  was  a  little  detailed.  The  use  of  a  graphics  tablet  attached 
to  an  IBM-AT,  however,  considerably  speeded  the  computation  of  both  the 
raingauge  polygons  and  soils  group  areas. 

6.2.1  Establishment  of  the  Data  Sets 


Two  data  sets  are  required  by  MILHY3;  'datal'  contains  the  program 
commands,  hydrological  commands  and  associated  data,  whilst  'data2’ 
contains  only  data  for  the  infiltration  algorithm.  The  rules  for  setting 
up  these  data  sets  and  examples  -f  them  are  given  in  Appendix  2. 

Datal 


Figure  6.3  illustrates  the  division  of  the  River  Fulda  catchment  into 
nine  subcatchments .  In  each  of  these  subcatchments,  a  runoff  hvdrograph 
must  be  developed,  and  in  all  except  the  headwater  subcatchments,  this 
must  be  added  to  the  flow  routed  through  the  subcatchment.  In  each 
routing  reach,  two  cross-sect  tons  are  developed,  one  at  either  end  of  the 
reach.  In  subcatchment  404,  where  the  River  Luder  joins  the  River  Fulda 
at  tts  inflow,  the  Kammerzell  cross-sect  ion  Is  used. 
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Table  6.3 


Soli  Group  Classification  for  the  Sub-Catchments 
on  the  Rivet  Fulda  Catchment 


Sub-catchment 

SC/SM 

UCSC  Soil 

Percentage 

ML 

Classification 

occurrence  of 

CH 

System 

group 

CL 

G 

401 

54.6 

11.6 

11.3 

10.5 

12.1 

402 

45.6 

10.3 

5.2 

27.7 

11  .4 

403 

25.0 

2.9 

4.0 

59.9 

7.9 

404 

36.6 

2.7 

15.2 

33.0 

12.0 

405 

65.8 

4.1 

4.7 

8.2 

17.3 

406 

50.1 

13.4 

9.8 

21.4 

5.4 

407 

46.4 

8.4 

25.2 

15.5 

4.6 

408 

41.0 

- 

15.2 

34.7 

9.2 

409 

86.5 

4.8 

8.6 

SC  -  Clayey  sands  or  clayey  gravelly  sands. 

SM  -  Silty  sand  or  silty  gravelly  sand. 

ML  -  Silts,  sandy  silts,  gravelly  silts. 

CH  -  Fat  clays 

CL  -  Lean  clays,  sandy  clays,  or  gravelly  clays 
G  -  Gravels  (grouped  together) 
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The  Curve  Numbers  for  the  generation  of  the  runoff  hydrograph  using  the 
SCS  method  were  identified  from  tables  in  the  Student  Handbook  on 
streamflow  forecasting  by  James  and  Stinson  (1981). 

Data2 


Each  soils  group  in  each  subcatchment  was  represented  by  a  single  soil 
column,  giving  a  total  of  42  columns  (see  Table  6.3).  The  runoff 
generated  by  these  columns  was  weighted  by  the  percentage  area  of  each 
subcatchment  that  a  soil  group  occupied.  For  the  six  storm  events 
identified  during  the  establishment  of  the  River  Fulda  catchment,  a 
common  theme  was  a  period  preceding  each  event  of  small  low  intensity 
showers.  This  enabled  the  fairly  safe  assumption  that  the  soils  were 
saturated  at  the  beginning  of  each  of  the  six  identified  events. 

Each  of  the  42  soil  columns  was  split  into  3  layers,  typical  of 
well-developed  soils,  and  a  total  of  10  cells  were  specified.  The  soils 
hydrological  parameters  were  identified  from  the  empirical  charts  and 
regression  equations  developed  by  Brakensiek  and  Rawls  (1983),  and 
reported  in  Anderson  and  Howes  (1984). 

To  check  that  these  data  sets  correctly  represented  the  River  Fulda 
catchment,  a  test  simulation  was  established  using  a  storm  occurring  in 
March  1986.  During  these  trials  two  problems  were  identified.  These 
were : 

i)  during  the  application  of  the  multiple  routing  technique,  it 

appeared  as  if  the  addition  subroutine,  'ADHYD',  that  sums  two 
hydrographs,  failed  to  be  operating  properly 

ii)  the  infiltration  algorithm  did  not  seem  to  be  generating  enough 

runoff  in  terms  of  the  total  volume  of  runoff  throughout  the  storm 

Investigation  into  the  first  of  these  problems  identified  a  small  bug  in 
the  original  coding,  that  enabled  one  of  the  time  intervals  of  the 
hydrographs  being  added  to  be  overwritten  if  the  time  intervals  of  the 
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two  hydrographs  were  different,  and  the  identity  number  of  the  sum  of  the 
two  hydrographs  was  the  same  as  one  of  the  hydrographs.  This  problem  was 
only  identified  because  of  the  complexity  in  the  structure  of  the  River 
Fulda  catchment,  and  the  large  number  of  additions  required.  It  is 
unlikely,  therefore,  that  any  previous  work  has  been  affected  by  this 
error.  This  bug  has  been  rectified  in  the  latest  version  of  MILHY3,  a 
copy  of  which  is  in  Appendix  3. 

The  second  problem  to  be  identified,  that  the  infiltration  algorithm 
failed  to  produce  enough  runoff,  was  firstly  blamed  on  the  low 
precipitation  intensity.  The  rainfall  for  each  subcatchment  was 
generated  by: 

1)  Identifying  the  daily  totals  for  each  of  the  rain  gauges  in  a 
Subcatchment . 

2)  These  totals  were  then  weighted  by  the  area  of  the  Theisson 
polygon  associated  with  each  raingauge. 

3)  The  daily  rainfall  for  each  subcatchment  was  then  divided 
into  three  time  intervals  of  eight  hours  each  according  to 
known  rainfall  figures  from  a  gauge  at  Fulda. 

4)  The  hourly  rainfall  was  then  computed,  assuming  a  minimum 
rainfall  intensity  of  lmm/hour,  distributing  backwards  from  the 
end  of  each  of  the  thrice  daily  intervals. 

5)  The  hourly  totals  were  then  cumulated. 

However,  this  method  and  assumptions  it  was  based  upon,  seemed 
reasonable,  as  the  volume  of  precipitation  over  the  River  Fulda  catchment 
was  conserved.  The  runoff  problem  seemed  to  be  a  soils  problem. 

As  noted  earlier,  the  soil  at  the  commencement  of  the  simulation  is 
assumed  to  be  saturated,  however  in  this  example  the  simulation  was 
started  a  more  realistic  three  days  prior  to  the  start  of  the  observed 
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outflow  hydrographs,  to  allow  some  drainage.  Analysis  of  the  results 
from  the  infiltration  algorithm  showed  that  most  of  the  precipitation  was 
draining  into  the  soil  and  very  little  was  generating  overland  flow. 

These  results  focused  attention  on  the  firstly  saturated  hydraulic 
conductivity  for  each  of  the  soil  classification  groups  and,  secondly, 
the  occurrence  of  each  of  these  soils  groups. 

A  second  trial  simulation  used  the  minimum  saturated  conductivity  values 
recorded  in  the  B-akensiek  and  Rawls  (1983)  tables  for  each  of  the  soils 
group.  This  did  not  increase  the  volume  of  runoff  generated ,. thereby 
identifying  the  soils  classif ication  maps  as  a  possible  source  of  error. 


6.2.2  Soils  classification  problems 

There  are  several  feasible  sources  of  error  in  the  generation  of  the 
proportion  of  a  subcatchment  that  a  soil  group  contributes.  These 
include : 

1.  Resolution.  This  includes  the  resolution  of  the  soil  surveyor's 
report,  and  the  interpretive  work  carried  out  in  the  establishment 
of  the  data  set.  The  resolution  of  the  original  survey  will 
depend  upon  the  purpose  to  which  the  map  is  aimed.  Beckett  and 
Webster  (1971)  point  out  that  there  is  little  practical  purpose  in 
having  a  resolution  size  less  than  the  minimum  land-use  management 
area,  usually  a  field. 

2.  Purity.  This  is  the  percentage  area  of  each  group  that  is 
occupied  by  that  group.  Beckett  and  Webster  (1971)  identify  the 
level  of  purity  many  of  the  soil  survey  organisations  attempt  to 
work  to,  and  this  includes  the  USDA  purity  level  of  80-90% 
(Simonsen,  1962),  and  the  US  Bureau  of  Reclamation  purity  level  of 
757.. 


Analysis  of  the  results  from  the  infiltration  algorithm  simulations 
showed  that  it  was  only  the  clay  groups  ( CH  and  CL,  see  Table  6.3)  that 
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had  a  low  enough  saturated  hydraulic  conductivity  to  generate  overland 
flow.  Could  we  legitimately  increase  the  percentage  occurrence  of  these 
two  groups,  in  line  with  the  purity  percentages  reported  above? 

Returning  to  the  soils  classification  maps,  closer  analysis  showed  that 
the  distribution  of  the  clay  grouped  was  heavily  biased  to  the  floodplain 
areas,  where  runoff  could  be  expected  to  join  the  channel  flow,  due  to 
the  high  water  table  and  low  slope  angles.  Without  entering  the  field  of 
"partially-contributing  areas",  and  models  such  as  VSAS-2  (Hewlett  and 
Troendle,  1975),  it  would  be  impossible  to  quantify  the  importance  of 
these  contributing  areas.  However,  it  would  suggest  that  in  future 
applications  the  percentage  area  of  soil  groups  with  low  hydraulic 
conductivities  should  be  estimated  more  accurately  at  the  expense  of 
accuracy  in  other  topographical  zones  of  the  catchment.  Having  completed 
the  measurement  of  each  soil  group's  extent  in  the  River  Fulda  catchment, 
it  was  decided  to  experiment  with  the  percentage  area  of  the  clay  groups 
in  the  three  headwateer  subcatchments. 

Trials  showed  that  increasing  the  occurrence  area  of  the  clay  groups  by 
35%  increased  the  runoff  volume  to  match  the  observed  hydrograph.  The 
consistency  of  the  level  of  increase  required  in  the  three  subcatchments 
confirmed  that  we  were  investigating  a  real  error  and  not  merely 
calibrating  the  model.  A  35%  increase  in  the  occurrence  area  of  the  clay 
groups  meant  that  at  least  15%  of  the  error  could  be  assigned  to  a 
resolution  problem,  accepting  that  the  purity  error  accounted  for  up  to 
20%  of  the  error.  Looking  at  the  soils  classification  maps  and  the 
complexity  of  soil  crops  in  the  valleys,  a  15%  error  seemed  very 
feasible.  In  the  simulations  reported  in  the  rest  of  Section  6, 
therefore,  the  increased  clay  percentage  occurrences  have  been  applied. 

6 . 3  Design  of  the  sensitivity  analysis 


The  main  objectives  of  the  sensitivity  analysis  are  to  investigate: 


1) 


The  sensitivity  of  the  hydrograph  to  variation  in  MILHY3's 
parameters . 
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2)  The  sensitivity  of  the  hydrograph  to  the  submodels  used  in  its 
development . 

3)  The  interaction  of  the  effects  of  the  turbulent  exchange  and 
multiple-routing  submodels. 

4)  The  relative  impact  of  the  inclusion  of  infiltration  algorithm  in 
catchments  where  subcatchments  are  utilized. 

5)  The  effects  of  scale. 

It  was  first  necessary,  therefore,  to  identify  all  the  parameters  and 
submodels  that  make  up  MILHY3.  Table  6.4  lists  all  the  variables  in 
three  groups,  spatial  variables,  temporal  variables  and  physically-based 
parameters.  The  spatial  and  temporal  variables  describe  the  resolution 
of  information  in  each  submodel  or  process  area,  whilst  the 
physically-based  parameters  describe  the  initial  conditions  and  geometry 
of  the  catchment.  Figure  6.7  describes  the  submodels  contained  in  MXLHY3 
and  shows  the  alternative  pathways  through  these  submodels. 

Table  6.5  contains  a  list  of  all  reported  sensitivity  analysis  carried 
out  on  all  versions  of  MILHY  and  HYMO.  The  table  also  gives  the  authors 
of  these  analysis  so  that  the  reader  may  refer  to  these  texts,  as  the 
design  of  the  sensitivity  analysis  reported  here  will  avoid  repetition  of 
this  work.  As  noted  earlier,  however,  the  simulations  reported  by 
Anderson  and  Howes  (1984,  1986),  and  in  Section  4  of  this  report,  are 
single  subcatchment  applications.  The  sensitivity  analysis  reported  here 
will  consider  the  relative  sensitivity  of  some  of  these  variables  in  a 
multi-subcatchraent  application. 

If  we  were  to  consider  varying  all  the  parameters  listed  in  Table  6.4  in 
a  statistically  meaningful  way,  then  a  conservative  estimate  of  the 
approximate  number  of  simulations  required  would  be  640,  taking  an 
impressive  6000  hours  of  CPU  on  the  SUN  workstation.  The  estimate 
ignores  the  ability  of  MILHY3  to  incorporate  stochastic  variability  in 
five  ^.f  the  so’1  parameters.  It  is  important,  therefore,  to  carefully 


Figure  6.7  :  MILHY3 
submodel  components 
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Table  6.4 

Variables  in  MILHY3 


Spatial  variability 

Number  of  subcatchments 

Number  of  rain  gauges 

Number  of  soil  columns 

Number  of  soil  cells 

Thickness  of  soil  cells 

Number  of  valley  sections  in  reach 

Number  of  segments  in  cross-section 

Rating  curve  stage  increment 

Temporal  variability 
Rainfall  time  increment 

Infiltration  simulation  iteration  period 
Routing  time  interval 

Physically-based  parameters 
Initial  soil  moisture  content 
Saturated  moisture  content 
Suction  moisture  curve 
Saturated  hydraulic  conductivity 
Surface  detention  capacity 
Maximum  evaporation 
Watershed  area 
Watershed  elevation 
Main  channel  length 
Cross-sectional  geometry 
Slope,  channel  and  floodplain 
Routing  length 

Mannings  'n',  channel  and  floodplain 
Reservoir  outflow  storage 

Soil,  crop,  conservation  and  gradient  factors 
Precipitation  -  storm  duration,  Intensity 


Table  6.5 


Sensitivity  Analysis  carried  out  on  MILHY 

Author  Parameter  or  Variable 

Williams  (1975)  Routing  length 

Routing  time  interval 

Detention  capacity 
Suction  moisture  curve 
Saturated  moisture  curve 
Saturated  hydraulic  conductivity 
Initial  moisture  content 

Anderson  &  Howes  (1986)  Watershed  area 

Watershed  elevation 
Main  channel  length 
Storm  intensity 
Cell  size 

Infiltration  simulation  increment 

Anderson  &  Singleton  (1987)  Precipitation,  spatial  variability 

Initial  results 
this  report: 

Section  4  Number  of  segments  in  cross-section 

Rating  cuive  stage  increment 
Cross-sectional  geometry 
Mannings  'n' 


Anderson  (1982) 

Anderson  &  Howes  (1984) 
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select  a  series  of  simulations  that  will  enable  the  sensitivity  of  MILHY3 
to  be  ascertained  with  a  fair  degree  of  certainty.  Interest  has  been 
directed  in  this  project  to  the  modelling  of  the  frictional  effects  of 
over-bank  flow  or  flow  in  two-stage  channels.  The  sensitivity  analysis 
will  reflect  this  concentration  and  therefore  focus  in  this  area.  The 
analysis  will  investigate  the  structure  of  MILHY3  in  terms  of  its 
submodel  components  and  explore  variability  in  the  physically-based 
parameters  in  the  downstream  conveyance  submodels.  These  parameters  will 
include  the  slope  of  the  channel  and  floodplain,  the  routing  length  and 
Manning  'n'  values  for  the  channel  and  floodplain. 

Having  directed  the  sensitivity  analysis  to  a  manageable  area  of 
investigation,  the  next  decision  was  to  investigate  alternative 
approaches  available  for  exploring  this  area. 

6.3.1  Alternative  methods  of  undertaking  a  sensitivity  analysis 
McCuen  (1973)  defines  sensitivity  as: 

"the  rate  of  change  in  one  factor  with  respect  to  change 
in  another  factor" 

McCuen  points  out  that  it  is  the  failure  of  researchers  to  appreciate  the 
generality  of  this  definition  that  has  limited  the  application  of  the 
sensitivity  analysis  tool  to  the  final  verification  of  models.  Several 
authors  have  identified  the  utility  of  the  sensitivity  analysts  in  all 
stages  of  the  development  of  a  model  (McCuen,  1973,  1976;  Miller  et  al. 
1976;  Hornberger  and  Spear,  1981).  This  is  why  an  initial  analysis  was 
incorporated  in  the  identification  of  most  important  processes  in 
two-stage  channels  (Section  3)  and  in  the  development  of  the  submodels  to 
model  these  processes  (Sections  \  and  5). 

Jones  (1982)  identified  two  approaches  to  undertaking  a  sensitivity 
analysis,  a  deterministic  and  a  stochastic  methodology.  A  deterministic 
methodology  Involves  making  small  changes  in  the  input  parameters  and 
investigating  changes  these  changes  provoke  on  the  model's  output.  A 
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stochastic  approach  involves  selecting  the  input  parameter  values  from  a 
probability  distribution  of  usually  either  the  accuracy  of  the  input 
values  or  the  error  bands  to  which  the  model  internally  operates.  A 
stochastic  analysis  can  usually  therefore  encompass  a  wider  range  of 
input  data  values  than  a  deterministic  analysis,  as  the  probability 
distributions  can  contain  all  physically  realistic  data  input  values.  In 
a  deterministic  analysis  the  operator  must  usually  either  select  input 
values  systematically,  or  use  his  intuitive  knowledge  of  the  behaviour  of 
the  input  parameters. 

Work  by  Anderson  and  Howes  (1984,  1986)  has  concentrated  on  a  stochastic 
analysis  of  the  soil  input  parameters  of  the  infiltration  algorithm.  The 
difficulties  of  measuring  these  parameters  in  the  field  mean  that  a 
relatively  large  error  band  can  be  associated  with  the  observed  field 
values.  A  stochastic  analysis  was  an  ideal  method  of  incorporating  these 
error  bands  in  an  analysis  of  the  effects  of  parameter  variability  on  the 
model's  output. 

In  the  sensitivity  analysis  reported  here,  we  ate  aiming  to  investigate 
the  behaviour  of  the  model's  output  with  respect  to  the  model  structure 
and  variability  in  certain  parameters.  As  the  model  structure  cannot  be 
described  by  a  probability  distribution,  this  necessitates  a 
deterministic  analysis. 

There  are  two  methods  of  computing  the  sensitivity  of  a  model  to  a 
parameter,  known  as  the  sensitivity  function,  in  a  deterministic 
approach.  These  are: 

1)  Differentiation:  the  model  described  as  a  function  is 
parametrically  differentiated  with  respect  to  each  parameter.  The 
mathematics  of  this  approach  has  yet  to  be  made  portable  enough  to 
enable  this  approach  to  be  widely  used. 

2)  Factor  perturbation:  each  parameter  is  incremented  and  the 
changes  in  the  output  quantified.  This  method  was  used  by  Smith 
(1976)  and  as  noted  earlier  has  extensive  computing  time 
requirements . 
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In  the  sensitivity  analysis  reported  here,  the  factor  perturbation 
technique  was  used  in  two  separate  approaches,  for  the  model  structure 
and  for  the  parameter  variability.  As  the  number  of  feasible  submodel 
combinations  was  relatively  small  (see  Figure  6.7),  it  was  possible  to 
establish  and  execute  a  complete  matrix  of  possible  permutations .  In 
investigating  the  sensitivity  of  MILHY3  to  variability  in  the  selected 
physically-based  parameters,  it  was  decided  to  attempt  to  utilize  u 
technique  previously  only  used  to  optimize  the  fit  of  parameters  in 
calibrated  models. 

Although  these  optimization  techniques  are  well-established,  it  seems 
that  they  have  not  been  used  as  an  alternative  to  the  sometimes  tedious 
development  of  factor  perturbation  matrices.  It  was  hoped  that,  if 
successful,  the  intermediate  values  of  the  optimization  scheme  would 
provide  a  good  indication  of  the  sensitivity  of  the  outflow  hydrograph, 
as  the  scheme  searched  to  find  the  'best-fit'  between  an  observed  and  a 
computed  hydrograph.  This  would  remove  the  necessity  to  initialise  a 
great  number  of  datasets  and  manually  search  through  the  results  sets.  A 
post-processor  could  search  through  the  iterations  of  optimization  scheme 
and  find  the  parameter  values  that  caused  the  greatest  or  smallest  impact 
on  the  computed  outflow  hydrograph.  Although  there  would  not  be  any 
direct  control  over  the  parameter  values  selected  by  the  optimization 
scheme,  in  a  factor  perturbation  analysis  the  selection  of  values  is 
usually  subjective  and  therefore  could  just  as  easily  overlook, 
significant  points.  However,  it  must  be  noted  that  this  investigation 
into  the  utility  of  optimization  schemes  is  rather  exploratory.  The 
feasibility  of  using  optimization  schemes  as  an  alternative  to 
traditional  factor  perturbation  techniques  will  depend  upon: 

1)  the  Initialization  period  to  set-up  the  scheme 

2)  Computer  demands 

1)  CPU 

ii)  disk  storage 
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6.3.2  The  optimization  approach 


As  noted  earlier,  the  technique  of  optimizing  the  fit  of  parameters  in 
hydrologic  models  using  a  sensitivity  analysis  for  the  purposes  of 
calibration,  is  well  established.  Applications  have  included  Armstrong 
et  al.  (1980)  and  Ibbit  and  O'Donnell  (1971).  McCuen  (1973)  identifies  a 
range  of  techniques  mostly  based  on  the  work  of  Cauchy  (1847),  who 
developed  the  method  of  converging  the  solution  utilizing  the  rate  of 
descent  or  gradient  of  an  objective  function  of  the  model's  output  in 
response  to  parameter  input  variability. 

A  range  of  optimization  techniques  for  minimising  and  maximising  a 
function  is  available  in  the  NAG  (Numerical  Algorithms  Group)  Library. 
Depending  on  the  level  of  sophistication  required  and  the  availability  of 
the  derivatives  of  the  function,  an  appropriate  routine  can  be  selected. 

A  simple  routine  was  selected  for  this  exploratory  investigation  (e04jaf) 
which  allows  the  user  to  select  the  upper  and  lower  boundaries  of  each 
variable  and  did  not  require  derivatives.  The  routine  wor’ ed  by 
developing  a  surface  of  values  for  a  function  (F)  that  describes  the 
difference  between  a  computed  value  and  an  observed  value.  The  routine 
then  searched  for  a  minimum  in  this  surface  by  selecting  parameter  values 
within  the  specified  boundaries.  A  prerequisite,  therefore,  of  this 
approach  is  that  a  function  can  be  computed  that  adequately  describes  the 
difference,  in  this  case,  between  an  observed  and  computed  hydrograpn. 

The  'least  squares'  approach  was  identified  as  being  a  function  already 
computed  by  MILHY3 ,  in  subroutine  'ERROR',  and  provided  a  simple  test  of 
the  fit  of  the  observed  hydrograph.  Figure  6.3  describes  how  the  MILHY3 
raodei,  the  function  and  optimization  routine,  e04jaf  fit  together 
schematically.  In  terras  of  the  computer  coding,  MLHY3  is  treated  as  a 
function  called  by  e04jaf,  which  i_  Itself  called  by  a  short  front 
programme  which  sets  up  the  boundary  conditions.  Once  the  routine  Is 
running,  it  Is  difficult  to  interrupt  as  all  the  commands  are  issued  by 
the  library  routine,  e)4jaf. 

As  this  Investigation  was  exploratory  In  nature  and  because  of  the 
concentration  of  the  analysts  on  the  downstream  conveyance  subroutines, 


the!  infiltration  algorithm  and  Curve  Number  routines  were  not  included  in 
the  optimization  scheme.  The  demands  of  the  processor  of  the  iterative 
nature  of  the  optimization  scheme,  and  storage  of  the  results  files  was 
foreseen  as  a  potential  problem  area,  without  exacerbation  from  the 
infiltration  algorithm. 

Setting  up  the  optimization  scheme,  shown  in  Figure  0.8,  proved  a 
reasonable  straightforward  task,  complicated  only  by  the  intermittent 
nature  of  the  'read'  statements.  The  'read'  statements  were  rewritten, 
so  that  all  commands  and  data  was  read  in  the  'front  end'  routine.  All 
'write'  statements  were  edited  out,  bar  the  warning  and  failure 
statements,  and  the  printing  of  the  outflow  hydrograph.  The  ability  of 
MILHY3  to  tolerate  any  set-up  structure  in  the  'datal'  dataset  was 
retained,  to  permit  the  use  of  the  multiple  routing  routine  which  is 
invoked  using  additional  commands  in  the  dataset. 

MILHY3  then  had  to  be  fronted  so  that  it  appeared  as  a  'function'  to  the 
optimization  routine  (e04jaf).  This  necessitated  the  addition  of  several 
COMMON  BLOCKS  to  ensure  all  the  data  was  correctly  passed  from  the 
initialization  (front-end)  routine.  Lastly,  all  the  parameters  had  to  be 
defined  as  being  'double-precision'  to  enable  them  to  be  correctly 
incremented  by  e04jaf.  It  took  about  two  weeks  to  set-up  and  test  the 
optimization  scheme,  a  similar  period  to  the  time  taken  to  set  up  a  whole 
series  of  datasets. 

To  test  that  the  optimization  scheme  was  working  properly  and  reaching  a 
minimum,  for  one  particular  application,  three  simulations  were 
undertaken.  In  each  of  these  simulations  the  initial  parameter  values 
were  changed  to  check  that  the  scheme  was  stopping  at  an  absolute 
minimum.  The  initial  conditions  specified  for  all  the  variables  were: 

1)  the  upper  boundary  limits 

2)  the  centre  of  the  parameter  limits 

3)  the  lower  boundary  limits 
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Figure  6.8  :  Conceptual  Optimization  Scheme 
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The  function  values  at  which  these  three  simulations  stopped  at  was 
however,  not  the  same.  Analysis  of  the  iterative  solution  showed  that 
each  of  the  simulations  had  got  stuck  in  local  minimums  close  to  the 
initial  conditions  in  each  case,  and  that  the  predicted  outflow 
hydrograph  and  parameter  values  at  each  of  the  local  minimums  were  widely 
different.  A  closer  examination  of  the  parameter  values  as  each 
simulation  converged  on  its  local  minimum  showed  that  the  parameters  were 

only  changing  from  iteration  to  iteration  by  a  very  small  amount, 

-4 

approximately  10  .  The  resolution  of  such  parameter  changes  was  too 

small  to  cause  any  change  in  the  predicted  outflow  hydrograph,  and  hence 
there  was  no  change  in  the  sum  of  squares  function  and  the  solution 
converged.  Three  problems  associated  with  the  resolution  of  parame.er 
variability  were  therefore  identified: 

1)  that  the  simulation  failed  to  converge  on  an  absolute  minimum 

2)  the  parameter  variability  increments  could  not  be  resolved 
with  the  accuracy  of  data  available  in  an  ungauged  catchment 

3)  the  parameter  variability  increment  caused  no  interpretable 
changes  to  the  predicted  hydrograph 

As  the  objective  of  this  investigation  was  to  pinpoint  parameter  changes 
that  did  or  did  not  cause  a  noticeable  effect  on  the  outflow  hydrograph, 
the  scheme  was  unsuccessful.  That  the  scheme  did  not  reach  an  absolute 
minimum  was  not  so  important  as  we  were  not  attempting  to  calibrate  the 
model.  The  reason  for  the  fail. re  to  reach  an  absolute  minimum  was 
important,  however,  as  it  showed  that  the  scheme  was  not  searching  within 
a  wide  enough  range  within  the  boundary  limits.  The  cause  of  this 
problem  was  simple;  it  was  the  size  of  the  parameter  changes  from 
iteration  to  iteration. 

A  solution  was  easily  come  by  -  a  different  optimization  scheme,  e04jbf, 
which  allows  the  user  to  select  the  resolution  of  parameter  changes. 

This  scheme  also  allowed  the  maximum  number  of  iterations  to  be  specified 
and  the  likely  size  of  the  sura  of  squares  at  the  absolute  minimum  to  be 
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estimated.  This  potentially  made  the  computing  demands  more 
controllable,  a  distinct  advantage,  the  main  disadvantage  being  that  the 
scheme  was  more  complex  to  set  up.  Having  already  made  the  program 
changes  necessary  in  MILHY3,  however,  the  bulk  of  the  work  was  already 
done  and  the  experience  of  setting  up  the  previous  scheme  meant  the 
changes  to  the  new  scheme  took  around  a  day  to  complete.  The  logic  of 
the  optimization  set-up  remained  unchanged  from  e04jaf,  shown  in  Figure 
6.8. 


Testing  the  new  optimization  scheme  from  a  variety  of  initial  conditions, 
resulted  in  solutions  with  a  discrepancy  well  within  the  level  of 
accuracy  to  which  the  parameters  could  be  estimated.  Initial 
investigation  into  the  iterative  results  proved  promising  enough  to 
encourage  us  to  pursue  the  application  of  this  optimization  scheme,  using 
data  from  the  River  Fulda  dataset. 

6.3.3  Application  of  the  factor  perturbation  technique 
to  the  Fulda  dataset 


From  the  analysis  of  the  possible  approaches  available  reported  above,  it 
can  be  seen  that  the  sensitivity  analysis  has  been  divided  into  two 
parts : 

1)  traditional  factor  perturbation  techniques  will' be  applied  to 
investigate  the  effects  of  the  model  structure 

2)  exploratory  optimization  techniques  will  be  used  to  investigate 
the  effects  of  variability  in  the  downstream  conveyance 
parameters . 

In  the  first  section  of  the  analysis,  two  events  will  be  explored,  an 
observed  out-of-bank  event  occurring  in  March  1986,  which  has  an 
occurrence  interval  of  1  in  10  years,  and  an  observed  In-bank  event  which 
occurred  in  June  1984,  These  two  events  will  be  simulated  for  the  whole 
of  the  River  Fulda  catchment  to  Rotenburg,  using  the  infiltration 
algorithm  and  the  curve  number  routine. 
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incorporating  turbulent  exchange  and  the  multiple  routing  routine  will 
also  be  used.  Throughout  these  simulations,  however,  the  only  parameters 
that  will  change  are  those  directly  related  to  the  storm  event.  In 
particular,  the  parameters  involved  in  the  downstream  conveyance 
submodels  will  remain  constant  throughout. 

In  the  second  exploratory  part  of  the  sensitivity  analysis,  the 
optimization  technique  will  be  used  to  test  the  sensitivity  of  the 
outflow  hydrograph  to  variability  in  five  parameters.  These  are: 

i)  floodplain  Mannings  'n' 

ii)  channel  Mannings  'n' 

iii)  floodplain  slope 

iv)  channel  slope 

v)  floodplain  routing  reach  length 

As  this  investigation  is  rather  exploratory,  a  single  storm  event  and  one 
river  reach  length  were  selected  for  exploration.  The  storm  event,  as 
above,  occurred  in  March  1986,  and  as  in  the  initial  analysis  reported  in 
Section  5,  the  reach  between  Bad  Hersfeld  and  Rotenburg  on  the  River 
Fulda  has  been  selected. 

The  results  from  the  investigation  into  the  impact  of  model  structure  are 
reported  in  Section  6.4,  whilst  the  optimization  results  are  reported  in 
Section  6.5. 


6 . 4  Sensitivity  of  the  outflow  hydrograph  to  model  structure 

The  objective  of  this  analysis  is  to  Investigate  the  sensitivity  of  the 
outflow  hydrograph  to  the  submodels  utilized  to  generate  it.  This 
analysis  will  also  investigate  the  effects  of  catchment  scale  on  the 
sensitivity  of  the  hydrograph. 

As  noted  earlier,  two  observed  storms  have  been  used,  the  first  occurred 
in  March  1986  and  has  a  recurrence  interval  of  1  in  10  years.  The  second 
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storm  occurred  in  June  1984  and  has  a  recurrence  interval  in  1  in  1.5 
years.  For  both  of  these  storms,  observed  hydrographs  were  available  for 
all  of  the  stations  shown  in  Figure  6.1,  except  for  the  outflow 
hydrograph  at  Rotenburg  for  the  storm  occurring  in  June  1984.  Both  of 
these  storms  are  out-of-bank  in  the  downstream  half  of  the  River  Fulda 
catchment . 

The  precipitation  data  is  used  by  the  curve  number  routine  in  MILKY,  and 
the  infiltration  algorithm  in  MILHY2  and  3  to  generate  Hortonian  runoff 
excess.  To  this  base flow  must  be  added.  For  the  purpose  of  this 
sensitivity  analysis  the  baseflow  levels  have  been  taken  from  the 
observed  hydrographs,  although  we  acknowledge  that  this  would  not  be 
feasible  in  an  ungauged  catchment.  There  is  potential,  therefore,  for 
the  development  of  a  submodel  to  generate  baseflow  conditions  based  on 
either  the  catchment  characteristics  or  channel  geometry  data  that  is 
already  derived. 

6.4.1  Storm  1  :  1  In  10  year  event 

Figure  6.9  shows  the  precipitation  pattern  at  Fulda  and  the  observed 
flood  hyetograph  at  Bad  Hersfeld  for  the  1  in  10  year  event.  The 
precipitation  patterns  temporal  distribution  for  all  the  subcatchments  is 
based  on  this  hyetograph,  and  the  magnitude  of  the  rainfall  in  each 
subcatchment  is  determined  from  the  daily  precipitation  records  of  the 
stations  shown  in  Figure  6.2.  The  minimum  rainfall  total  for  the  event 
occurred  in  subcatchraent  408  where  45  mm  fell,  the  maximum  occurring  in 
subcatchment  403  where  75mm  fell.  The  observed  hydrograph  shown  in 
Figure  6.9  illustrates  that  the  discharge  peak  occurred  approximately  24 
hours  after  the  rainfall  peak.  The  time  to  peak  of  the  observed 

hydrograph  is  30  hours  from  the  simulation  start  time  and  the  peak 
3  -1 

discharge  is  426  m  s 

The  predicted  outflow  hydrographs  at  Bad  Hersfeld  using  the  curve  number 
and  the  infiltration  algorithm  are  recorded  in  Tables  6.6  and  6.7 
respectively.  These  tables  summarize  the  characteristics  of  the 
hydrographs,  noting  the  peak  discharge,  time  to  peak  and  equivalent 
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Table  6.6 


Computation 

Method 

Peak  discharge 

3-1 
ra  s 

Time  to  peak 

hours 

Runoff  depth 

ra 

IT-1 

MR-0 

265 

17.5 

0.03 

IT=2 

MR-0 

272 

17.5 

0.03 

IT- 3 

MR-0 

265 

17.5 

0.03 

IT-4 

MR =0 

249 

17.5 

0.03 

IT-1 

MR-1 

312 

16.5 

0.03 

IT-2 

MR-1 

328 

16.5 

0.03 

IT- 3 

MR-1 

323 

16.5 

0.03 

IT-4 

MR-1 

281 

17.0 

0.03 

IT  = 

turbulent  exchange 

routine 

MR  = 

multiple  routing  - 

invoked  =  1 

- 

not  invoked  =  0 

IT-1 

vertical  interface 

,  zero  shear 

IT-2 

vertical  interface 

,  apparent  shear  stress  ratio  =  1 

IT- 3 

diagonal  interface 

,  zero  shear 

IT- 4 

diagonal  interface 

,  apparent  shear  stress  ratio  =  1 
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Table  6.7 

Storm  1  :  1  In  10  year  event 


Predicted  Outflow  at  Bad  Hersfeld  utilizing  the 
Infiltration  Algorithm 


Computation  Method 

Peak  discharge 

mV1 

Time  to  peak 

hours 

Runoff  depth 

m 

IT-1 

MR-0 

312 

18.5 

0.03 

IT-2 

MR-0 

321 

18.0 

0.03 

IT- 3 

MR-0 

310 

18.0 

0.03 

IT-4 

MR-0 

290 

18.0 

0.03 

IT-1 

MR-1 

364 

16.5 

0.04 

IT-2 

MR-1 

383 

16.5 

0.04 

IT-3 

MR-1 

372 

16.5 

0.04 

IT-4 

MR-1 

332 

17.5 

0.04 

IT  =  turbulent  exchange  routine 

MR  =  multiple  routing  -  invoked 

=  1 

not  invoked  =  0 
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runoff  depth,  for  all  the  turbulent  exchange  and  multiple  routing  routine 
combinations . 

Figure  6.10  illustrates  the  impact  of  the  computation  method  on  the 
predicted  outflow  hydrograph  at  Bad  Hersfeld.  MILHY  invokes  the  curve 
number  routine,  whilst  MILHY2  and  3  utilize  the  infiltration  algorithm. 
MILHY  and  MILRY2  use  the  turbulent  exchange  method  2,  and  single  routing 
reaches,  whilst  MILHY3  uses  exchange  method  3  and  the  multiple  routing 
reach  routine. 

Comparison  of  Tables  6.6  and  6.7  shows  that  both  the  curve  number  and 
infiltration  algorithm  produce  runoff  depths  comparable  with  the  observed 
data  for  the  drainage  area  up  to  Bad  Hersfeld.  Analysis  of  the  runoff 
volumes  generated  at  the  intermediate  stations  (see  Figure  6.1)  shows 
that  this  is  true  for  all  the  subcatchments.  Tables  6.8  and  6.9 
illustrate  the  runoff  volumes  for  the  station  at  Hermannspiegal ,  where 
the  observed  runoff  volume  is  0.028  metres. 

Comparison  of  the  curve  number  routine  and  infiltration  algorithms 
prediction  of  the  peak  discharghe  rate  at  Bad  Hersfeld  and  Hermannspiegal 
(Figures  6.10  and  6.11)  shows  that  in  both  cases  the  infiltration 
algorithm  produces  higher  peak  discharges.  This  behaviour  has  been  noted 
previously  by  Anderson  (1982)  and  Anderson  and  Howes  (1984),  who  showed 
that  this  behaviour  occurred  during  high  and  low  intensity  storms.  In 
this  analysis  it  is  worth  noting  that  this  behaviour  is  still  visible 
after  the  hydrographs  have  been  routed  through  up  to  four  subcatchraents. 
This  re-emphasises  earlier  work  by  Anderson  and  Howes  (1984,  1986)  that 
illustrated  the  importance  of  the  shape  of  the  runoff  hydrograph  in 
determining  the  outflow  hydrograph. 

Analysis  of  Tables  6.6  to  6.9  shows  that  the  peak  discharge  is  the 
parameter  most  sensitive  to  the  downstream  computation  method.  The 
impact  of  the  turbulent  exchange  and  the  multiple  routing  routines  is 
related  to  the  depth  of  flow  on  the  floodplains.  The  outflow  hydrograph 
at  Hermannspiegal  (shown  in  Figure  6.11)  is  routed  from  Marbach.  The 
distribution  of  the  inflow  hydrograph  at  Marbach  means,  however,  that 
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Table  6.8 

Storm  1  :  1  In  10  year  event 


Predicted  Outflow  at  Hermannsplegal  utilizing  the 
Curve  Number  Routine 


Computation  Method  Peak  discharge 

3  -1 
m  s 

Time  to  peak 

hours 

Runoff  depth 

m 

IT-1 

MR-0 

70 

13.5 

0.02 

IT- 2 

MR-0 

70 

15.0 

0.02 

IT- 3 

MR-0 

70 

15.0 

0.02 

IT-4 

MR-0 

69 

15.0 

0.02 

IT-1 

MR-1 

70 

15.0 

0.02 

IT-2 

MR-1 

71 

15.0 

0.02 

IT- 3 

MR-1 

69 

15.0 

0.02 

IT-4 

MR-  l 

68 

15.0 

0.02 

IT  = 

turbulent  exchange 

routine 

MR  = 

multiple  routing  - 

invoked  - 

1 

- 

not  invoked  = 

0 
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Table  6.9 


Storm  1  :  1  In  10  year  event 

Predicted  Outflow  at  Hermannspiegal  utilizing  the 
Infiltration  Algorithm 


Computation  Method 

Peak  discharge 

3  -1 
m  s 

Time  to  peak 

hours 

Runoff  depth 

m 

IT-1 

MR-0 

89 

16.0 

0.03 

IT-2 

MR-0 

90 

16.0 

0.03 

IT- 3 

MR-0 

90 

16.0 

0.03 

IT- 4 

MR-0 

89 

16.0 

0.03 

IT-1 

MR-1 

90 

16.0 

0.03 

IT-2 

11 

X 

90 

16.0 

0.03 

IT- 3 

73 

II 

89 

16.0 

0.03 

IT-4 

MR-1 

87 

16.0 

0.03 

IT  =  turbulent  exchange  routine 
MR  =  multipLe  routing  -  invoked  =  1 
-  not  invoked  =  0 


L 


JL 


145 


only  three  data  points  at  the  peak  of  the  hydrograph  are  assigned  to  the 
floodplains.  The  impact  on  the  turbulent  exchange  and  multiple  routing 
routines  is  therefore  minimal.  At  Bad  Hersfeld,  the  impact  of  the  new 
routines  Is  more  pronounced.  Tables  6.6  and  6.7  show  that  the  turbulent 
exchange  techniques  that  assume  zero  shear,  methods  1  and  3,  produce  very 
similar  results  and  the  methods  that  assume  an  apparent  shear  stress 
ratio  of  1  produce  smaller  discharge  predictions.  These  results  are  as 
expected,  as  the  methods  that  assume  zero  shear  do  not  incorporate  the 
shear  surfaces  between  segments  in  out-of-bank  flow  and  therefore  the 
type  of  division  has  nominal  effect.  The  methods  that  incorporate  the 
shear  stresses  surfaces  would  be  expected  to  have  lower  discharge 
predictions,  as  the  larger  wetted  perimeter  reduces  the  hydraulic  radius 
and  hence  the  discharge  using  the  Manning  equation.  The  discharge 
discrepancy  between  turbulent  exchange  remains  relatively  constant, 
irrespective  of  the  absolute  discharge. 

All  of  the  simulation  methods  predicted  the  time  of  the  hydrograph  peak 
too  early.  This  can  be  related  to  the  rainfall  data  which  was  derived 
from  records  taken  at  eight  hour  intervals. 

6.4.2  Storm  3:  1  in  1.5  year  event 

Figure  6.12  shows  the  hydrograph  for  subcatchment  402  and  the  observed 
discharge  hydrograph  at  Bad  Hersfeld  for  the  1  in  1.5  year  event.  The 
precipitation  data  for  each  subcatchment  was  derived  from  hourly  data 
available  for  the  station  at  Bad  Hersfeld.  The  spatial  variability 
was  generated  from  the  daily  records,  and  showed  a  minimum  total 
precipitation  in  subctachment  403  of  58  mm,  and  a  maximum  in  subcatchment 
406  with  71  mm.  In  contrast  to  the  1  in  10  year  event,  this  event  is 
characterized  by  a  double  peak  in  the  rainfall  event,  the  effects  of 
which  can  be  seen  in  the  observed  hydrograph.  In  both  of  these  peaks, 
the  discharge  peak  occurs  approximately  30  hours  after  the  precipitation 
peak. 

The  results  from  the  curve  number  and  infiltration  algorithm  simulations 
are  given  in  Tables  6.10  and  6.1 L  for  Bad  Hersfeld,  and  6.12  and  6.13  for 
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Table  6.10 

Storm  3  :  1  In  1.5  year  event 

Predicted  Outflow  at  Bad  Hersfeld  utilizing  the 
Curve  Number  Routine 


Computation 

Method 

Peak  discharge 

mV1 

1  2 

Time  to  peak 

hours 

1  2 

Runoff  depth 

m 

IT-1 

MR-0 

238 

355 

19.5 

55.5 

0.04 

IT-2 

MR-0 

242 

359 

19.5 

55.0 

0.04 

IT-3 

MR-0 

235 

361 

19.5 

55.0 

0.04 

IT-4 

MR-0 

224 

353 

19.5 

55.5 

0.04 

IT-1 

MR=1 

262 

364 

19.0 

54.5 

0.04 

IT- 2 

MR-1 

273 

383 

19.0 

54.0 

0.04 

IT- 3 

MR-1 

273 

377 

19.0 

54.5 

0.04 

IT-4 

MR-1 

224 

352 

19.5 

55.5 

0.04 

IT  =  turbulent  exchange  routine 
MR  =  multiple  routing  -  Invoked  =  1 
-  not  invoked  =  0 
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Table  6.11 


Storm  3  :  1  In  1.5  year  event 

Predicted  Outflow  at  Bad  Hersfeld  utilizing  the 
Infiltration  Algorithm 


Computation  Method 

Peak  discharge 

3  -1 
m  s 

1  2 

Time  to  peak 

hours 

1  2 

Runoff  depth 

m 

IT-1 

MR-0 

236 

238 

16 

55.0 

0.04 

IT-2 

MR-0 

241 

232 

16 

54.5 

0.04 

IT-3 

MR-0 

238 

235 

16 

54.5 

0.04 

IT-4 

MR-0 

229 

235 

16 

54.5 

0.04 

IT-1 

MR-1 

240 

255 

16 

54.0 

0.04 

IT-2 

MR-1 

249 

269 

16 

54.0 

0.04 

IT-3 

MR-1 

253 

266 

16 

54.5 

0.04 

IT-4 

MR-1 

226 

234 

16 

54.0 

0.04 

IT  =  turbulent  exchange  routine 
MR  =  multiple  routing  -  invoked  =  1 
-  not  Invoked  =  0 
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Figure  6.14 
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Table  6.12 

Storm  3  :  1  In  1.5  year  event 

Predicted  Outflow  at  Hermannsplegal  utilizing  the 
Curve  Number  Routine 


Computation 

Method 

Peak  discharge 

■V* 

1  2 

Time  to  peak 

hours 

1  2 

Runoff  depth 

m 

IT-1 

MR-0 

81 

138 

18.0 

52.0 

0.05 

IT-2 

MR-0 

81 

139 

18.0 

52.0 

0.05 

IT-3 

MR-0 

81 

138 

18.0 

52.0 

0.05 

IT-4 

MR-0 

81 

135 

18.0 

52.5 

0.05 

IT-1 

MR-1 

81 

138 

18.0 

52.0 

0.05 

IT-2 

MR-1 

81 

140 

18.0 

52.0 

0.05 

IT- 3 

MR-1 

81 

132 

18.0 

52.0 

0.05 

IT-4 

MR-1 

81 

134 

18.0 

52.0 

0.05 

IT  =  turbulent  exchange  routine 
MR  =  multiple  routing  -  invoked  =  1 
-  not  invoked  =  0 
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Table  6.13 


Storm  3  :  1  In  1.5  year  event 

Predicted  Outflow  at  Bad  Hersfeld  utilizing  the 
Infiltration  Algorithm 


Computation 

Method 

Peak  discharge 

3  -1 
m  s 

Time  to 

hours 

peak 

Runoff  depth 

m 

1 

2 

1 

2 

IT-1 

MR-0 

67 

88 

13 

52 

0.05 

IT-2 

MR-0 

67 

88 

13 

52 

0.05 

IT- 3 

MR-0 

67 

89 

13 

51 

0.05 

IT-4 

MR-0 

66 

88 

13 

52 

0.05 

IT-1 

MR-1 

66 

88 

13 

52 

0.05 

IT-2 

MR-1 

66 

89 

13 

52 

0.05 

IT- 3 

MR-1 

67 

87 

13 

52 

0.05 

IT-4 

MR-1 

66 

84 

13 

52 

0.05 

IT  =  turbulent  exchange  routine 
MR  =  multiple  routing  -  Invoked  =  1 
-  not  Invoked  =  0 
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Hermannspiegal .  These  tables  record  the  peak,  discharge  and  time  to  peak 
of  both  of  the  peaks  and  the  runoff  depth  of  the  simulation.  The  results 
are  summarized  in  Figures  6.13  for  Bad  Hersfeld  and  Figure  6.14  for 
Hermannspiegal . 

Comparing  the  curve  number  and  infiltration  algorithm  simulations  (MILHY 
and  MILHY2) ,  in  all  the  subcatchments  the  curve  number  routine  generates 
far  too  large  a  discharge  in  the  second  peak.  In  the  infiltration 
algorithm,  the  effects  of  the  second  stage  of  the  precipitation  on  the 
discharge  are  negated  by  the  preceding  dry  period  during  which  drainage 
occurs.  In  the  infiltration  algorithm,  therefore,  some  of  the 
precipitation  is  absorbed  by  the  soil  columns  before  saturated  conditions 
reoccur  and  overland  flow  is  predicted.  This  more  complex  storm 
therefore  shows  the  superior  predictive  capability  of  MILHY2  over  MILHY, 
in  multiple  and  single  subcatchment  applications. 

Analysis  of  the  impact  of  the  turbulent  exchange  and  multiple  routing 
routines,  shows  as  in  the  1  in  10  year  the  peak  discharge  is  the 
hydrograph  most  sensitive  to  any  change  in  the  computation  method.  The 
impact  at  Hermannspiegal  is  negligible  because  the  upstream  hydrograph  is 
in-bank.  At  Bad  Hersfeld  the  greatest  impact  is  generated  by  the 
multiple  routing  routine. 

6.4.3  Comparison  of  the  two  storms 

Figures  15  and  16  illustrate  the  impact  of  the  multiple  routing  and 
turbulent  exchange  routines  on  the  predicted  outflow  hydrograph  at  Bad 
Hersfeld  for  both  of  the  storms.  It  is  interesting  to  compare  the 
effects  of  these  routines  individually  and  their  impact  when  operating 
together  which  is  illustrated  by  the  MILHY3  simulations  in  Figures  6.10 
and  6.13.  The  turbulent  exchange  method  (IT)  used  throughout  this 
analysis  is  method  3  (diagonal  Interface  with  zero  shear  stress),  as  the 
results  recorded  in  Tables  6.7  and  6.11  show  this  method  seems  to  have 
the  greatest  impact  on  the  peak  discharge. 

Figure  6.16  shows  that  the  impact  of  the  turbulent  exchange  method  3,  has 
no  noticeable  impact  on  the  predicted  hydrographs  of  either  storm,  when 
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Figure  6.16 
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applied  without,  the  multiple  routing  routine.  Figure  6.15  shows, 
however,  that  the  multiple  routing  routine  has  a  significant  effect 
on  the  hydrograph,  and  that  the  routine's  impact  varies  between  the  two 
storms.  In  Storm  1  the  main  peak  of  the  hydrograph  is  reduced  whilst  the 
minor  peaks  on  the  recession  limb  are  accentuated.  In  Storm  2,  by 
contrast,  the  recession  from  the  first  peak  is  steepened  and  the  second 
peak  is  significantly  accentuated. 

Taking  each  of  the  storms  in  turn,  in  Storm  1  it  is  important  to 
appreciate  that  floodplain  flow  only  occurs  from  Unter-Schwarz  to  Bad 
Hersfeld  on  the  River  Fulda.  The  inflow  hydrographs  at  Unter-Schwarz  for 
the  simulation  with  multiple  routing  routine  are  identical  to  that 
without  the  routine,  and  that  the  inflow  from  the  Haune  is  channel  flow 
only  and  therefore  can  be  ignored.  Comparison,  therefore,  concentrates 
on  the  travel  time  tables  for  the  Unter-Schwarz  to  Bad  Hersfeld  reach  for 

the  two  simulations.  At  the  main  peak  of  the  hydrograph  of  around  350 
3  -1 

ms  ,  in  the  application  without  multiple  routing  the  time  taken  for  the 
peak  to  travel  the  length  of  the  reach  is  13  hours.  In  the  multiple 
routing  application,  however,  45%  of  this  peak  is  apportioned  to  the 
floodplain  where  the  travel  time  is  approximately  19  hours.  The 
remaining  55%  is  assigned  to  the  main  channel  where  the  travel  time  is 
only  s  hoiir<s .  This  difference  in  travel  times  means  that  the  peak  of  the 
hydrograph  is  flattened  out.  In  the  later  minor  peaks,  the  effect  of  the 
division  of  floodplain  and  channel  flows  is  rather  different.  Looking  at 
two  points  on  the  inflow  hydrograph,  the  travel  time  without  multiple 
routing  is  11.4  hours,  with  multiple  routing  the  floodplain  travel  time 
is  55  hours,  and  the  main  channel  travel  time  is  6  hours.  However,  as 
only  4%  of  the  flow  is  assigned  to  the  floodplain,  with  multiple  routing 
the  flow  arrives  earlier  and  the  peaks  are  more  accentuated.  Storm  1 
shows,  therefore,  that  the  effects  of  the  multiple  routing  routine  on  the 
outflow  hydrograph  is  determined  by  the  percentage  of  flow  that  is 
assigned  to  the  floodplain. 


> 


; 


Storm  3  confirms  this  conclusion  as  where  15%  of  the  flow  is  assigned  to 
the  floodplain,  the  multiple  routing  prediction  is  more  attenuated. 

Where  floodplain  flows  account  for  10%  or  less  of  the  total  discharge, 
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the  hydrographs  less  attenuated,  as  the  multiple  routing  channel  time  is 
significantly  lower  than  the  joint  channel/ floodplain  travel  time. 

The  predicted  hydrograph  at  Bad  Hersfeld  for  the  application  of  both  the 
multiple  routing  and  turbulent  exchange  routines  (MILHY3)  are  shown  in 
Figures  6.10  and  6.13  for  Storms  1  and  3  respectively.  Comparison  with 
Figures  6.15  and  6.16  shows  that  the  hydrograph  for  Storm  1  is 
significantly  different  from  those  produced  by  the  multiple  routing  and 
turbulent  exchange  routines  alone.  The  Storm  3  hydrograph  exactly 
matches  the  hydrograph  from  the  multiple  routing  routine.  It  would  seem, 
therefore,  that  the  effects  of  applying  both  routines  varies  according  to 
the  storm. 

In  Storm  1,  the  MILHY3  prediction  seems  particularly  strange  as  when 
compared  to  the  MILHY2  solution,  MILHY3  predicts  the  main  peak  as  being 
earlier  and  attenuates  it  less.  This  contrasts  with  Figure  6.15,  where 
the  multiple  routing  routine  increases  the  attenuation  of  the  peak. 
Analysis  of  the  rating  curves  and  travel  time  tables  generated  by  MILHY3, 
and  those  from  the  simulation  shown  in  Figure  6.i5,  showed  that  it  was 
the  travel  times  that  control  the  attenuation  of  the  hydrograph.  When 
the  routines  are  applied  together,  the  travel  times  are  reduced  and  more 
flow  is  rssigned  to  the  floodplain.  The  turbulent  exchange  generated 
small  changes  in  the  rating  curve  and  changes  In  the  travel  time  that 
have  no  impact  on  the  hydrograph  when  applied  on  its  own,  hence  Figure 
6.16.  When  applied  with  the  multiple  routing  routine,  these  small 

changes  become  significant.  For  example,  at  t he  fifth  hour  of  the 

3  - 1 

simulation,  for  IT=2,  MR=l,  18%  (39  ms  )  of  the  total  discharge  was  in 
the  left  floodplain,  this  water  has  a  travel  time  of  70  hours.  In 
contrast,  when  IT=3,  MR=1  (MILHY3)  26%  (58  m^s  *)  of  the  discharge  was  in 
the  left  floodplain,  the  travel  time  was  60  hours. 

In  Storm  3,  however,  there  are  no  noticeable  differences  between  the 
multiple  routing  and  multiple  routing  with  turbulent  exchange  solutions. 
Although  there  are  differences  in  the  travel  time  tables  between  these 
two  techniques,  these  differences  do  not  become  significant  as  a  much 
greater  proportion  of  the  hvdrograph  is  out-of-bank. 
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6.4.4  Conclusions 

From  the  analysis  of  the  March  1986  and  June  1984  storms  it  is  possible 

to  conclude  that:- 

1)  The  infiltration  algorithm  still  makes  a  considerable  improvement 
in  the  predictive  capability  of  MILHY  when  several  subcatchments 
are  applied. 

2)  The  multiple  routing  routine  has  a  significant  impact  on  the 
predicted  hydrograph.  When  floodplain  inundation  accounts  for 
approximately  15 Z  or  more  of  the  total  discharge,  then  the 
multiple  routing  routine  increases  the  attenuation  of  the 
floodwave.  When  floodplain  flows  account  for  less  than  10%  of  the 
total  discharge,  the  multiple  routing  routine  decreases  the 
attenuation  of  the  floodwave. 

3)  The  turbulent  exchange  routine  makes  no  significant  impact  on  the 
predicted  discharge  hydrograph. 

4)  When  the  multiple  touting  and  turbulent  exchange  routines  are 
applied  together,  then  changes  invoked  by  the  turbulent  exchange 
routine  becomes  significant.  When  15“  or  mor  of  the  discharge  is 
assigned  to  the  floodplain,  the  floodwave  attenuation  is  reduced. 

5)  For  cases  where  the  floodplain  accounts  for  15%  or  more  of  t lie 
discharge,  then  the  joint  application  of  the  multiple  routing  and 
turbulent  exchange  routines  improves  the  predictive  capability  of 
MILHY . 

6 . 5  Opt  imiz.it  ion  results 

\s  noted  earlier,  this  investigation  into  the  utility  of  optimization 

techniques  as  sensitivity  tools  is  rather  exploratory,  and  therefore  i 

single  reach  ind  storm  are  used.  The  reach  use.!  is  between  tile  gauging 
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stations  at  Bad  Hersfeld  and  Rotenburg  on  the  River  Fulda,  and  has  been 
described  in  Section  4.4.1.  An  observed  hydrograph  of  the  1  in  10  year 
event  occurring  in  March  1986  was  specified  as  the  upstream  inflow,  and 
the  computed  outflow,  assuming  r.o  inflow  from  the  intervening  drainage 
area,  was  compared  with  the  observed  hydrograph  at  Rotenburg. 

Two  methods  of  computing  F,  the  function  that  describes  the  difference 
between  the  computed  and  observed  hydrographs,  were  used,  to  see  if  this 
had  any  impact  on  the  range  of  parameter  values  the  optimization 
technique  would  select.  The  techniques  were: 

i)  ordinary  sum  of  squares  6.1 

n 

0F^  *  ^  ^  CVmr'CVrf1c)'' 

i-l 


ii)  absolute  error  divided  by  variance 


where  qm 
qc 
qm 


0 FI  -r  - i*i 


2  (y^-vT0,i 


^  C^rnc-^ no )a 


CM 


observed  peak,  discharge 
computed  peak  discharge 
-  mean  observed  discharge 


6.: 


Both  of  these  methods  of  analysis  are  incorporated  in  Lhe  subroutine 
"ERROR",  which  remains  unchanged  from  the  version  incorporated  in  MILHY2. 

Tlie  upper  and  lower  boundary  limits  for  the  five  parameters  identified 
earlier  are  reported  in  Table  6.14,  which  also  includes  the  initial 
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Table  6.14 

Initial  Conditions,  Boundary  Conditions  and  Variable  Increments 
for  the  Optimization  Simulations 


Initial  Boundary  Limits  Variable 

Conditions  Upper  Lower  Increments 


Floodplain 

Mannings  'n' 

0.05 

0.16 

0.025 

0.01 

Channel 

Mannings  'n' 

0.035 

0.1 

0.025 

0.01 

Floodplain 

Slope 

0.0006 

0.001 

0.0001 

0.0001 

Channel 

Slope 

0.0007 

0.001 

0.0001 

0.0001 

Floodplain 

21951 

23750 

16860 

1525 

routing  reach 
length  (m) 


161 


values  for  the  values  at  the  start  of  each  simulation,  and  the  variable 
increment  intervals. 

The  number  of  iterations,  MAXCAL,  was  set  randomly  at  500.  All  of  the 
simulations  reached  this  maximum  without  satisfying  the  optimization 
routines  rules  for  a  minimum,  and  consequently  failed.  Investigation 
into  the  results  files,  however,  showed  that  the  routines  rules  for  a 
minimum  were  very  stringent  and  for  the  level  of  accuracy  required  in 
this  analysis,  an  effective  minimum  was  reached  at  around  400  iterations. 

As  the  computing  demands  of  this  approach  were  foreseen  as  being  a  major 
drawback  in  the  utility  of  the  technique,  a  close  eye  was  kept  on  the  CPU 
and  size  of  files  produced.  Output  at  each  iteration  was  restricted  to 
the  parameter  values  and  the  function  (F)  value,  and  the  post-processor 
added  the  computed  outflow  hydrograph  approximately  every  tenth 
iteration.  CPU  demands  for  a  single  reach  varied  from  400  to  800  seconds 
for  500  iterations.  Trial  simulations  for  the  whole  Fulda  catchment, 
utilizing  the  curve  number  routine  to  generate  runoff,  took  up  to  6000 
seconds  of  CPU.  For  applications  where  CPU  demands  are  small, 
therefore,  the  optimization  technique  is  an  efficient  method  of 
undertaking  multiple  simulations.  With  the  utilization  of  the 
infiltration  algorithm,  the  CPU  demands  for  the  Fulda  catchment  would 
become  excessive,  as  each  simulation  iteration  would  take  9  hours  of  CPU. 

The  results  files  proved  to  be  more  of  a  problem  as  the  files  were  very 
large.  Sizes  ranged  from  12  kbytes  to  2.9  mega  bytes,  files  too  large  to 
edit.  Clearly  the  post-processor  needs  to  be  more  selective  in  the 
iteration  results  It  saves.  Because  of  '"he  size  of  the  results  files  the 
results  presented  here  represent  only  a  small  selection  of  the  most 
important  points. 

The  presentation  of  the  results  has  been  structured  in  order  to  answer 
several  questions.  These  are: 

i)  to  which  of  the  five  parameters  is  the  hydrograph  most  sensitive? 

i.i)  does  the  computation  method  affect  this  sensitivity? 
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iii)  do  the  parameters  interact  to  increase  or  decrease  the  sensitivity 
of  the  hydrograph? 

In  much  of  the  analysis  presented,  the  relative  error  and  absolute  error 
are  used  to  compare  the  differences  between  the  computation  methods.  The 
relative  error  is  dimensionless  and  is  defined  as: 

RE  =  x  -  x .  6.3 

c  l 

x . 

l 

where  xc  -  computed  value  of  x 

x^  -  value  of  x  under  initial  conditions 

6.5.1  Sensitivity  to  parameter  variability 

Tables  6.15  to  6.19  show  the  absolute  (AE)  and  relative  error  (RE)  for 
the  peak  discharge,  time  to  peak  and  sum  of  squares  for  five  of  the 
computation  methods  and  all  five  of  the  variables.  The  errors  shown  are 
computed  for  one  increment  step  above  and  below  the  initial  conditions. 

Analysis  of  these  tables  show  the  asymmetrical  sensitivity  of  the  three 
indicators,  peak  discharge,  time  to  peak  and  sura  of  squares,  around  the 
initial  conditions.  For  example,  Table  6.15  shows  that  the  sensitivity 
of  the  peak  discharge  to  variation  in  the  floodplain  Mannings  'n'  value 
is  markedly  different  for  values  greater  than  the  initial  conditions  than 
values  less  than  the  initial  conditions.  This  asymmetrical  effect  is 
particularly  noticeable  for  the  variation  in  Mannings  'n',  both  in  the 
floodplain  and  channel  (see  Tables  6.15  and  6.16),  suggesting  that  the 
sensitivity  of  the  hydrograph  to  variation  in  'n'  is  not  linear. 

Tables  6.15  to  6.19  show  the  sensitivity  of  the  hvdrograph  to  a  one 
incremental  step  in  the  mid-point  between  the  upper  and  lower  boundary 
limits  for  all  five  variables.  Table  6.20  shows  an  example  of  the 
relative  errors  generated  from  a  one  increment  step  increase  in  each  of 
the  parameters  at  the  boundaries  and  compares  these  with  the  mid-limit 


Time  to  peak 


Increment 


Computation 

method 


Peak  discharge 


$?- 


RE 


AE 

hours 


Sum  of  squares 
0F2 

RE  AE  RE 


i 


IT=2 

MR=0 

-125 

0.01 

0 

0.00 

62941 

0.00 

IT=3 

MR=0 

-115 

0.13 

0 

0.08 

62951 

0.31 

+1 

IT=I 

MR-1 

-107 

0.01 

+9 

0.00 

124796 

0.12 

IT-2 

MR-1 

-127 

0.10 

+9 

0.07 

105431 

0.12 

IT-3 

MR-1 

-  88 

0.03 

0 

0.00 

86893 

0.08 

IT-2 

MR-0 

-121 

0 

63177 

IT-3 

MR-0 

-  72 

-3 

91568 

0 

IT-1 

MR-1 

-110 

+9 

111038 

IT-2 

MR-1 

-  95 

+6 

94106 

IT-3 

MR-1 

-  79 

-3 

94782 

IT- 2 

-123 

0.01 

0 

0.00 

68801 

0.09 

IT- 3 

MR-0 

-  71 

0.00 

-3 

0.00 

92675 

0.01 

-1 

IT-1 

MR-1 

-  71 

0.13 

+6 

0.06 

92118 

0.17 

IT-2 

MR-1 

-102 

0.02 

+3 

0.07 

70169 

0.25 

IT- 3 

MR-1 

-  70 

0.03 

-3 

0.00 

106950 

0.13 
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Table  6.16 


Errors  from  one  Increment  step  variation  in  channel  Mannings  'n' 


Increment 

Computation 

method 

Peak  discharge 

$E_!  RE 

ms 

Time  to 

AE 

hours 

peak 

RE 

Sum  of  squares 
0F2 

AE  RE 

IT=2 

MR-0 

-138 

0.06 

0 

0.00 

68270 

0.08 

IT-3 

MR-0 

-  91 

0.06 

-3 

0.00 

74276 

0.19 

+  1 

IT-1 

MR-1 

-143 

0.11 

+9 

0.00 

85134 

0.23 

IT-2 

MR-1 

-103 

0.03 

+9 

0.07 

107817 

0.15 

IT- 3 

MR-1 

-  82 

0.01 

-3 

0.00 

91518 

0.03 

0 

IT-2 
IT- 3 
IT-1 
IT- 2 
IT- 3 

MR-0 

MR-Q 

MR=l 

MR-1 

MR-1 
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-  72 
-no 

-  95 

-  79 

0 

-3 

+9 

+6 

-3 
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111038 

94106 

94782 

IT- 2 

MR-0 

-105 

0.06 

0 

0.00 

68694 

0.09 

IT- 3 

MR-0 

-  96 

0.07 

-3 

0.00 

73639 

0.19 

+  1 

IT-1 

MR-1 

-  96 

0.05 

+3 

0.13 

70847 

0.36 

IT-2 

MR-1 

-  58 

0.12 

+3 

0.07 

68635 

0.27 

IT- 3 

MR-1 

-  73 

0.02 

-3 

0.00 

107051 

0.13 

Errors  from  one  Increment  step  variation  In  floodplain  slope 


Increment 

Computation 

method 

Peak  discharge 

^E  RE 

m  s 

Time  to 

AE 

hours 

peak 

RE 

Sum  of  squares 
0F2 

AE  RE 

IT=2 

MR=0 

-120 

0.00 

0.0 

0.00 

63421 

0.00 

IT=3 

MR=0 

-  72 
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-3 

0.00 
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+  L 

IT=1 
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-105 

0.02 

+9 
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0.02 

IT-2 

MR=  1 

-  89 
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+6 

0.00 
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IT- 3 
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-  75 
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-3 
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0 
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MR-1 

MR-1 
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-3 

+9 

+6 

-3 
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94783 

IT- 2 

MR-0 

-123 

0.01 

0.0 

0.00 
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0.00 

IT- 3 

MR-0 

-112 

0.12 

0.0 

0.08 

63756 

0.30 

-1 

IT-1 

MR-1 

-  99 

0.04 

+9 

0.00 

124144 

0.12 

IT- 2 

MR-1 

-110 

0.05 

+3 

0.07 

104471 

0.11 

IT- 3 

MR-1 

-  83 

0.01 

-3 

0.00 

90447 

0.05 

Errors  from  one  increment  variation  la  channel  slope 


Increment 

Computation 

method 

Peak  discharge 

RE 

m  s 

Time  to 

AE 

hours 
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RE  ' 

Sum  of  squa.as 
0F2 

AE  RE 

IT=2 

MR-0 
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0.00 
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MR=0 
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0.10 

0  .no 

0.08 

66138 

0.28 

+  1 

IT-1 

MR=  1 

-112 

0.01 

+9 

0.00 

107855 

0.03 

IT- 2 

MR=1 

-100 

0.02 

+6 

0.00 

90519 

0.04 

IT=3 

MR-1 

-  78 

0.00 

-3 

0.00 

96431 

0.02 

IT-2 

MR-0 

-121 

0.0 

63177 

IT- 3 

MR-0 

-  72 

-3 
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0 

IT-1 
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+9 
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IT-2 

MR-1 

-  95 

+6 

94106 

IT- 3 

MR-1 

-  79 

-3 

94732 

IT- 2 

MR-0 

-126 

0.02 

IT- 3 

MR-0 

-  78 

0.02 

IT-1 

MR-1 

-103 

0.02 

IT-2 

MR-1 

-121 

0.08 

IT-3 

MR-1 

-  80 

0.00 

0.00 

0.00 

63223 

0.00 

-3 

0.00 

84485 

0.08 

+9 

0.00 

110171 

0.01 

+9 

0.07 

108635 

0.15 

-3 

0.00 

93286 

0.02 

-1 
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Table  6.19 


Errot <s  frog  one  Increment  step  variation  in  floodplain  routing  reach  length 


Relative 

error  in  peak  discharge 

from  one  increment 

variation 

at  the  boundaries 

and  mid-way 

Variable 

Upper  Boundary 

Relative  Errors 
Mid-range 

Lower  Boundary 

Floodplain 
Mannings  'n' 

0.053 

0.10 

0.01 

Channel 

Mannings  'n' 

0.003 

0.025 

0 

Floodplain 

slope 

0.083 

0.019 

0.058 

Channel 

slope 

0.083 

0.014 

0.03 

Routing  reach 
length 


0.029 


0.014 


0.011 
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values.  This  particular  example  compares  the  relative  errors  in  the  peak 
discharge  for  one  computation  method  and  is  typical  of  the  tables  derived 
for  other  computation  methods. 

Table  6.20  shows  that  at  both  the  upper  and  lower  boundaries  the  outflow 
hydrograph  is  most  sensitive  to  slope.  In  the  mid-ranges,  the  hydrograph 
is  most  sensitive  to  Mannings  'n'.  For  MILHY  applications  significant 
inundation  therefore  is  important  to  define  Mannings  ' n'  as  accurately  as 
possible;  this  is  especially  true  in  the  floodplain.  Table  6.20  also 
shows  that  the  outflow  hydrograph  is  not  sensitive  to  relatively  small 
changes  in  the  floodplain  routing  length,  except  when  slopes  are  steep  (1 
x  10~3). 

6.5.2  Sensitivity  variations  associated  with  the  computation  method 

The  computation  method  of  the  optimization  scheme  incorporates  two 
sources  of  variation:- 

1)  the  structure  of  MILHY3,  specifically  in  this  case  which  of  the 
turbulent  exchange  routines  has  been  used  and  whether  the  multiple 
routing  routine  has  been  invoked 

2)  the  factor  (F)  used  to  quantify  the  differences  between  the 
observed  and  computed  hydrographs 

Section  6. A  investigated  the  sensitivity  of  the  outflow  hydrograph  to 
variations  in  MILHY3's  structure  for  two  storms  in  the  River  Fulda 
catchment.  Here  the  impact  of  the  model  structure  on  the  sensitivity  of 
the  hydrograph  to  parameter  variability  is  investigated.  Tables  6.15  to 
6.19  compare  the  relative  errors  for  three  measures  of  hydrograph  fit, 
for  a  range  of  model  structures.  If  the  sensitivity  of  the  model  to 
variation  in  the  five  physically  based  parameters  were  not  affected  by 
the  model  structure  then  we  would  expect  the  relative  errors  for  all  the 
computation  methods  to  be  the  same.  The  fact  that  there  are  variations 
suggests  that  certain  computational  techniques  increase  the  sensitivity 
of  the  model  to  parameter  change.  This  problem  is  particularly 


noticeable  in  the  sensitivity  to  variation  in  floodplain  Mannings  ' n'  and 
in  the  channel  Mannings  'n'  with  the  introduction  of  the  multiple  routing 
routine . 

The  second  question  raised  in  this  section  is:  does  the  function  (F) 
utilized  to  describe  the  fit  of  the  predicted  affect  the  utility  of  using 
optimization  techniques  as  part  of  a  sensitivity  analysis?  Analysis  of 
the  results  showed  that  the  exact  function  value  did  not  influence  the 
routine's  selection  of  parameter  values;  only  the  relative  function 
value  between  simulations  was  used.  The  solutions  from  both  functions 
converged  on  minimums  for  which  the  five  parameter  values  were  very 
close.  We  accept,  however,  that  this  may  not  be  the  case  if  a  radically 
different  function  were  applied. 

6.5.3  Conclusions 


The  optimization  esults  have  shown  that:- 

1)  Optimization  schemes  provide  a  viable  alternative  to  traditional 
factor  perturbation  sensitivity  analyses  provided  that: 

(a)  the  CPU  demands  of  the  model  can  be  met 

(b)  a  satisfactory  function  can  be  found  to  describe  the  fit  of 
the  hydrograph 

2)  The  sensitivity  of  MILHY3  in  two-stage  applications  is  dominated 
by  :- 

(a)  slope  when  slopes  are  ^  1  x  If)  ^ 

(b)  floodpLain  Mannings  'n'  when  slopes^  1  x  10  4 

3)  Slopes  need  to  be  defined  to  an  accuracy  of  1  x  10  ^ 

i)  Mannings  'n'  values  need  to  be  defined  to  an  accuracy  of  at  least 

1  x  10-2 
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5)  The  short-circuiting  of  floodplain  flow  is  only  significant  if  the 

-3 

floodplain  slope  is  4  1  x  10  and  the  floodplain  routing  length 
is  at  least  10%  shorter  than  the  main  channel  routing  length. 


Conclusions 


In  answer  to  the  specific  questions  posed  at  the  start  of  Section  6, 

it  is  possible  to  conclude  from  this  analysis  that:- 

1)  The  sensitivity  of  the  predicted  outflow  hydrograph  to  variation 
in  its  parameters  and  process  submodels  is  dependent  on  the 
magnitude  and  duration  of  the  out-of-bank  event. 

2)  All  out-of-bank  events  are  sensitive  to  the  inclusion  of  the 
multiple  routing  routine.  Larger  out-of-bank  events  are  sensitive 
to  the  joint  application  of  the  turbulent  exchange  and  multiple 
routing  routine,  where  MILHY3  significantly  improves  the  accuracy 
of  the  predicted  hydrograph  over  MILHY2. 

_3 

3)  MIIHY  is  sensitive  to  variations  in  slope  J  1  x  10  and 

_  9 

floodplain  Mannings  ' n'  variations  of  ^  1  x  10  '. 

!*)  Variability  in  the  predicted  hydrograph  can  be  attributed  to  the 
slope  parameter  when  slopes  are  ^  0.001,  and  Mannings  'n'  when 
slopes  <  0.001  and ^  0.0001. 

3)  The  introduction  of  the  infiltration  algorithm  is  still  the  most 

important  improvement  in  the  predictive  capability  of  MILHY. 
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FLOODPLAIN  MODELLING 


7.1  Introduction 

The  preceding  chapters  have  been  concerned  with  the  establishment  and 
implementation  of  MILHY3.  As  we  have  discussed,  such  a  scheme  would 
appear  to  have  significant  advantages  over  the  earlier  versions  of  MIL'riY . 
Notwithstanding  this,  however,  the  MILHY  schemes  only  provide  for  the 
estimation  of  stage  at  subcatchment  or  routing  reach  outflow  points. 

Considerable  benefit  would  accrue  if,  firstly,  more  within-reach 
information  were  available  in  a  predictive  context,  and,  secondly,  if  the 
method  to  generate  this  was  eventually  compatible  with  the  MILHY3 
capability  outlined  above. 


One  of  the  objectives  of  the  current  research  was,  therefore,  to 
establish  whether  an  existing  finite  element  method  ( FEM) ,  in  the  form  of 
RMA-2,  could  be  implemented  at  a  sufficiently  large  reach  length  (13- 30 
km).  As  we  will  outline  below,  this  scheme  is  capable  of  generating 
detailed  out  of  bank,  within  reach  predictions  of  stage  and  hence 
inundated  area.  Wo  sought,  therefore,  to  attempt  the  implementation  of 
the  scheme  on  the  River  Fulda  and  to  undertake  initial  validation  tests. 

7 . 2  Existing  Applications  )f  Finite  Element  Methods  for 
River  Reach  Studies 


Two  dimensional  horizontal  finite  element  numeri  :al  models  have  been 
applied  to  certain  clisses  of  river  channel  problems.  Applications  ha. 
included  detailed  analyses  of  flow  patterns  near  structures  such  is 
bridges  (Tseng,  19  7";  -ling  ami  Norton,  ,  river  conf luence  stu.il  .• 
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alternatives  to  physical  models  and  the  scale  of  interest  has  been  small, 
e.g.  reaches  of  river  a  few  river  widths  long.  An  exception  has  been 
some  estuary  studies  that  were  of  large  scale  (e.g.  tens  of  miles);  some 
of  which  utilized  a  "hybrid"  (numerical  plus  physical)  modelling 
technique  (McAnally  et  al.,  1984a  &  1984K'.  Moreover,  in  a  rev'ew  of  the 
application  of  finite  element  methods  tu  river  channels,  Samuels  reports 
that  in  only  two  studies  (Lee,  1980  and  Zielke  &  Urban,  1981)  is  the 
river  channel  resolved  separately  from  the  floodplain. 

7 .3  Large  Scale  Floodplain  Modelling  with  Finite  Element  Methods 

From  the  above  review  of  model  applications ,  the  lack  of  attention  to 
large  scale  floodplain  modeling  is  evident.  This  is  the  primary  focus  of 
this  chapter  and  is  drawn  from  two  separate  research  needs: 

( i)  to  establish  the  accuracy  of  finite  element  modelling  at  larger 
river  reach  scales  (of  the  order  of  16-32  km)  and  appropriate 
operational  rules  for  such  applications 

(ii)  to  explore  the  feasibility  of  using  finite  element  modelling  (at 
this  larger  scale)  to  effectively  extend  the  record  of  a  gauged 
catchment  to  include  extreme  events  which  may  not  be  available. 
These  large  events  may  be  needed  for  the  purposes  of  assessing  the 
accuracy  of  simpler  hydrological ly  based  forecasting  models.  The 
proviso  here  is  that  the  finite  element  method  is  considered  to  be 
capable  of  providing  icourate  predictions  ((i)  above)  that  cun  be 
equated  with  known  'gauged'  data  where  extreme  event  data  is 
unava liable . 

7.4  Model  Selection 

lue  numerical  mode!  known  is  R'fA-2  (\inc  in.!  '.'ortm,  197.9)  was  select*! 

:  >r  use  in  this  study  ml  ill  the  runs  ne  essiry  for  the  project  were 
eider*,  iken  by  hr  1)  '•!  dee  it  the  Myd r  ) in 1  -  Tag ineer in-;  renter,  Davis  , 
California.  This  m>!ei  solves  the  lepth  into  ;r  ited  Reynolds  equations 
:  ir  t  v-  i-d  inn  ns  i  '  n  i '.  '  r-'e-  s  a-  •  ,  ,  flow  in  the  ,r  i  t  i pi  lie  us  in  f 
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using  the  finite  element  technique.  It  can  be  applied  to  both  steady  and 
unsteady  flows.  The  finite  element  formulation  of  RMA-2  allows  boundary 
roughness  and  solution  resolution  to  vary  spatially  to  accurately  reflect 
topography.  It  also  provides  for  a  wide  variety  of  boundary  conditions, 
including  stage  hydrographs,  discharge  hydrographs,  and  rating  curves. 
Recent  research  has  produced  schemes  that  allow  elements  to  be  either  wet 
or  dry;  that  is,  to  have  a  horizontally  moving  flow  boundary. 

RMA-2  has  proved  a  useful  tool  for  solving  a  wide  variety  of  complex  flow 
problems  where  the  traditional  one-dimensional  assumption  that  the  flow 
path  and  distribution  can  be  determined  a  priori  is  questionable.  Most 
applications,  however,  have  been  to  either  large  estuaries  (McAnally,  et 
al,,  1984a;  McAnally  et  al.,  1984b,  MacArthur  et  al.,  1987)  or  to  the 
prediction  of  local  flow  patterns  -ear  structures  (Gee,  1985). 

The  ability  of  RMA-2  to  allow  dry  areas  within  the  solution  domain  during 
the  simulation  of  an  unsteady  flow  event  such  as  a  flood  wave  led  us  to 
select  it  for  testing  on  a  floodplain  problem  where  flow  is  initially 
within  the  channel,  spreads  into  the  overbank  areas  as  the  flood  arrives, 
and  returns  to  the  channel  as  the  flood  recedes.  The  two-dimensional 
solution  relieves  the  engineer  of  having  to  construct  cross  sections  that 
are  perpendicular  to  the  flow  for  all  flows,  as  is  required  in  a 
one-dimensional  analysis. 

The  version  of  RMA-2  (version  4)  used  in  this  study  contains  a  new 
approach  to  the  we t t ing/ dry ing  problem.  Previously,  an  element 
instantaneously  became  dry  once  the  depth  at  any  node  in  that  element 
became  zero  or  negative,  and  similarly  with  wetting.  This  resulted  in 
relatively  large  changes  in  cross-sectional  area  as  overbank  elements 
were  brought  into  the  solution.  The  new  approach  is  based  upon  the 
concept  of  "marsh"  elements  that  gradually  wet  and  dry.  This  is 
accompl ished  through  a  pseudo-porositv  that  operates  .an  the  flow  eirrving 
capacity  of  an  element  as  the  depth  changes  (King,  1987).  The 
application  described  herein  is  the  first  practical  application  of  the 
car  ala  element  option. 
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7 . 5  Research  Design 


In  exploring  Che  capability  of  finite  element  models  (such  as  RMA-2)  t 
operate  on  large  scale  floodplain  applications ,  there  are  two  principal 
approaches  that  could  be  taken  to  the  overall  research  design.  Firstly, 
a  single  reach  study  could  be  undertaken  utilizing  accurate  field  and 
recorded  data,  or,  secondly,  hypothetical  reaches  could  be  established 
and  estimated  parameter  values  usede  It  was  felt  desirable  in  this  case 
to  select  the  former  approach  because  we  wished  to  eventually  explore  the 
relative  effect  of  field  parameter  values  on  the  model  result,  in 
relation  to  the  geometrical  definition  of  the  application  as  well  as  the 
model  structure.  This  opportunity  would  have  been  denied  if  hypothetical 
reaches  were  used.  We  do,  however,  acknowledge  that  since  we  only 
utilize  the  model  on  one  study  reach,  our  conclusions,  while  hopefully 
appropriate  to  other  floodplains  with  similar  characteristics,  must  be 
interpreted  within  this  framework.  As  will  be  evident  from  the  following 
sections,  the  model  set-up  and  subsequent  computation  time  requirements 
prevented  a  greater  number  of  study  reaches  from  being  established  under 
this  research  project. 

The  study  reach  selected  for  the  RMA-2  application  was  that  of  the  Bad 
Hersfeld-Rotenburg  section  of  the  River  Fulda  in  West  Germany.  This  is  a 
24  km  reach  that  is  described  in  the  following  section,  and  is  of  a 
topographic  type  that  is  typical  not  only  of  the  region,  but  many  other 
river  systems  in  Europe. 

With  any  field  application  data  limitations  may  subsequently  become 
evident.  In  particular,  in  that  we  were  seeking  to  examine  the  potential 
of  finite  element  schemes  to  model  large-scale  floodplain  inundation,  we 
were  particularly  interested  in  high  magnitude  extreme  events.  While  1 
in  10  year  events  may  be  available  in  terms  of  observed  hydrographs,  1  in 
LOO  year,  or  larger,  events  may  not  be.  These  flows  are  of  interest  in 
the  current  context  because  of  the  associated  extreme  inundation  areas. 

It  is  likely,  therefore,  as  in  this  application,  that  high  magnitude 
events  need  to  be  generated  from  smaller  events  by  application  of 
accepted  flood  frequency  techniques.  This  in  no  way  represents 
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a  criticism  of  the  overall  approach  we  have  adopted,  in  that  for 
hypothetical  reach  studies,  totally  hypothetical  events  would  be  used 
with  'unknown'  recurrence  intervals.  In  the  context  of  data  acquisition, 
available  hydrographs  and  related  gauging  information  was  obtained  from 
the  Fulda  River  Authority. 

Four  simulations  were  undertaken  with  RMA-2  using  two  overbank 
roughnesses  for  each  of  two  flood  hydrographs.  To  develop  the  larger  of 
the  two  hydrographs  an  observed  1  in  10  year  flood  was  scaled  up  to  a 
1  in  100  year  flood.  Both  floods  were  modelled  using  field  estimated 
values  of  Manning's  n.  To  enable  an  evaluation  of  model  sensitivity, 
prediction  susceptibility  to  parameter  change  and  possible  error  in  field 
estimation  of  roughness,  the  two  overbank  roughness  conditions  were 
considered  to  be  essential  to  the  research  design. 

7 . 6  Study  Reach 


Having  established  that  the  study  was  to  be  based  upon  a  natural  river 
reach,  the  River  Fulda  between  Bad  Hersfeld  and  Rotenburg  was  selected 
(figure  5.1).  The  reach  is  some  15  miles  (24  km)  long  with  a  slope  of 
0.0008  (a  50  ft.,  or  15  m. ,  drop  in  elevation).  The  floodplain  is 
typically  0.6  miles  (1  km)  wide.  From  figure  2,  it  can  be  seen  that  the 
floodplain  itself  has  a  very  shallow  slope  orthogonal  to  the  river 
(typically  0.0001),  and  is  bounded  by  steep  hills,  often  forested.  The 
floodplain  land  use  is  typically  grazed  pasture.  Field  estimates  of 
Manning's  n  were  undertaken  throughout  the  reach  using  the  photographic 
definition  of  roughness  type  identified  in  Chow  (1959).  Figure  6.6 
illustrates  channel  and  floodplain  conditions  that  dominate  the  majaority 
of  the  reach  and  illustrates  an  important  category  of  flood  inundation 
problem  in  comparatively  sinuous  rivers,  especially  where,  as  in  this 
case,  both  settlements  and  transportt ation  routes  are  located  within  the 
floodplain  area.  Roughness  was  assessed  as  0.045  for  the  floodplain  and 
0.035  for  the  channel.  There  were,  however,  a  number  of  sections  within 
the  reach  where  somewhat  rougher  floodplain  conditions  were  evident,  as 
shown  in  figure  4.9.  These  latter  floodplain  areas  were  assessed  as 
having  a  Manning's  a  of  1.07. 
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Channel  cross-sectional  data,  rating  curve  and  hydrographs  are  available 
for  the  Bad  Hersfeld  and  Rotenburg  sections  and  were  used  in  the  RMA-2 
application  detailed  below.  At  bankfull,  the  channel  is  about  13  ft. 
deep  and  98  ft.  wide  (4m  by  30m)  at  Bad  Hersfeld  and  is  18  ft.  deep  and 
164  ft.  wide  (5.5m  by  50m)  at  Rotenburg. 

Available  maps  of  inundation  show  the  whole  floodplain  to  have  been 
inundated  in  1946.  However,  no  record  of  the  flood  hydrograph  for  this 
event  was  found. 


7.7  System  Schematizat ion 


RMA-2  utilizes  a  finite  element  network  composed  of  both  triangular  and 
quadrilateral  elements.  Ground  elevations  are  defined  at  the  corners  of 
the  elements  nd  assumed  to  vary  linearly  between  corner  nodes.  In  this 
study,  the  channel  was  represented  by  a  strip  of  two  elements  wide 
(figure  7.1)  producing  a  triangular  cross  section.  Overbank  areas  were 
represented  by  much  larger  elements.  The  lateral  extent  of  the  network 
was  determined  by  a  bluff  line,  beyond  which  none  of  the  simulated  flood 
events  would  extend. 

Manning's  roughness  coefficients  may  be  varied  spatially  in  RMA-2 
applications.  For  this  study,  two  roughnesses  were  used  in  each 
simulation,  one  for  the  channel  elements  and  another  for  the  overbank 
elements.  The  resulting  finite  element  network  was  composed  of  860 
elements  and  2660  nodes.  The  ratio  of  maximum  to  minimum  element  sices 
was  about  200  to  1.  This  variability  in  resolution  demonstrates  the 
flexibility  of  the  finite  element  method  for  use  in  large  scale 
floodplain  modelling. 

The  times  simulated  were  21  hours  for  Storm  1  and  34  hours  for  Storm  2. 
These  times  covered  the  rising  limb  of  the  hydrographs  and  a  portion  of 
the  recession  Limb.  AIL  results  reported  herein  were  accomplished  using 
a  0.5  hr.  time  step.  No  over- ittenuat  ion  due  to  this  relatively  large 
time  step  was  observed.  One  simulation  of  Storm  1  was  performed  using 
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an  0.25  hr.  time  step,  yielding  results  the  same  as  those  with  the  0.5 
hr.  time  step.  Further  investigations  need  to  be  done  regarding  the 
selection  criteria  and  sensitivity  of  results  to  various  time  steps  for 
use  of  this  implicit  finite  element  model  for  floodplain  inundation 
studies . 

7 .8  Computational  Aspects 


Although  this  is  not  a  very  computationally  intensive  problem  for  the 
simulation  of  steady  flow  conditions,  the  dynamic  simulations  performed 
(consisting  of  40  to  60  time  steps)  utilized  significant  computational 
resources.  Each  simulation  took  several  hours  of  central  processing  time 
on  a  Harris  100  super  microcomputer  at  H.E.C.  Although  contemporary  desk 
top  computers  equal  or  exceed  the  processing  speed  of  this  computer,  the 
results  indicate  that  engineers  contemplating  two-dimensional  floodplain 
modelling  on  this  scale  for  dynamic  flow  events  should  carefully  plan 
their  studies  to  minimize  the  number  of  situations  to  be  modelled  and 
utilize  steady  flow  simulations  wherever  possible.  It  should  be  pointed 
out  here  that  traditional  floodplain  Inundation  mapping  is  accomplished 
using  serai-empirical  discharge  routing  with  steady,  one-dimensional 
computation  of  water  surface  elevation  to  define  inundated  areas. 

7.9  Results 


As  discussed  above,  two  storms  were  applied  to  the  study  reach.  Storml 
was  an  event  commencing  31  March  1986,  which  produced  a  1  in  10  year 
flood  event  at  the  downstream  station  (Rotenburg).  Continuously  recorded 
stage  hydrographs  were  available  for  this  storm  at  both  Bad  Hersfeld  and 
Rotenburg.  Using  RjMA-2  in  association  with  the  mesh  illustrated  in 
figure  7.1,  Storm  1  was  simulated  using  values  of  Manning's  n  for  the 
channel  of  0.035  and  for  the  floodplaim  of  0.045.  The  results  of  this 
simulation  are  illustrated  in  figures  7.2a  and  7.2b.  From  these  results, 
It  is  evident  that  there  is  under-prediction  of  the  stage  at  the  upstrern 
location  (figure  7.2a)  and  over-prediction  of  stage  by  about  1  ft.  (0.3:n) 
at  the  peak  stage  at  the  downstream  station  (figure  7.2b).  These  results 
indicated  that  the  channel  plus  floodplain  conveyance  was  too  efficient. 

A  second  simulation  was  performed  based  upon  the  floodplain  condition 
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shown  in  figure  4.9  (overbank  n  of  0.07)  and  the  results  are  shown  on 
figure  7.3. 

Storm  2  was  a  generated  flood  event  using  the  pattern  of  Storm  1  ,  scaled 
up  to  a  1  in  100  year  peak  flow  using  available  flood  frequency  data  at 
Bad  Hersfeld  and  Rotenburg.  The  results  of  utilizing  the  0.045  and  0.035 
Manning's  n  coefficients  are  shown  in  figure  7.4.  These  results  show  a 
similar  trend  to  those  depicted  in  figure  7.2  (Storm  1)  in  terms  of 
predicted  values.  However,  in  terms  of  absolute  stage  prediction,  the 
results  for  Storm  2  show  a  greater  correspondence  with  the  'observed' 
stage.  This  result  reinforces  the  observation  we  made  above  relating  to 
floodplain  roughness  because  in  Storm  2,  a  significantly  greater  area  of 
the  floodplain  is  inundated.  At  peak  stage,  the  over-prediction  at 
Rotenburg  is  0.6  ft.  (0.2m),  with  the  arrival  of  the  peak  being  about  2 
hours  early  (see  figure  7.2b).  This  timing  error  corresponds 
approximately  with  the  resolution  at  which  it  is  possible  to  assess  the 
stage  from  the  gauge  recordings. 

Following  the  scheme  adopted  for  Storm  1,  a  second  simulation  of  Storm  2 
was  undertaken.  Figure  7.3  shows  the  simulations  corresponding  to  0.07 
and  0.035  roughness  values  for  the  floodplain  and  channel  respectively 
for  Storm  2  (see  figure  4.9).  Figure  7.5a  indicates  a  particularly  good 
correspondence  with  the  'observed'  stage  compare  to  figure  7.6a 
(floodplain  roughness  of  0.045).  In  addition,  an  improvement  in  the  peak 
stage  timing  (figure  7.3b)  was  observed. 

Intermediate  stage  hydrographs  were  also  examined  for  locations  with  n 
the  study  reach.  These  showed  complete  consistency  with  both  the 
upstream  and  downstream  -esults  reinforcing  the  suggestion  that  ttie  model 
may  indeed  be  used  to  estimate  the  extent  of  floodplain  inundation 
throughout  an  entire  reach  at  the  scale  used  in  this  study. 

The  two-dimensional  solution  obtained  from  RMA-2  produces  velocity 
vectors  in  addition  to  stage  at  every  computational  node.  Indeed  most 
applications  of  two-dimensional  clow  models  have  focused  on  velocity  for 
purposes  of  constituent  transport  or  design.  In  the  context  of  large 
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floodplain  modelling  velocities  are  important  for  both  definition  of 
inundated  area,  and  definition  of  flood  hazard.  Figure  7.6  shows  a 
typical  velocity  vector  plot  computed  for  the  peak  discharge  of  Storm  1 
(refer  to  figure  7.3).  For  this  1  in  10  year  event,  the  predictions 
shown  in  figure  7.6  imply  substantial  overbank  flow  throughout  the  entire 
reach.  Whilst  no  field  inundation  data  are  available  for  Storms  2  or  2 
for  validation  purposes,  this  prediction  is  consistent  with  the  limited 
published  field  evidence  relating  to  a  1946  flood  that  was  larger  than 
Storm  1.  This  fact,  taken  together  with  what  appears  to  be  a 
satisfactory  reproduction  of  the  stage  hydrographs  at  both  upstream  and 
downstream  locations  for  the  observed  Storm  1 ,  suggests  that  it  is 
reasonable  to  conclude  from  the  available  evidence  thag  the  predicted 
inundated  overbank  areas  (figure  7.6)  are  realistic. 

7.10  Discussion 

Application  of  RMA-2 ,  version  4,  to  the  River  Fulda  has  enabled  a 
preliminary  assessment  to  be  made  of  the  applicability  of  finite  element 
numerical  models  to  large  scale  floodplain  applications.  The  initial 
results  presented  in  this  paper  indicate  that  RMA-2  may  successfully  be 
used  for  estimating  the  depth  and  lateral  extent  of  inundation  at  this 
scale.  Also,  the  distribution  of  flow  velocities  across  the  floodplain 
is  available  from  the  simulations;  there  were  no  data,  however,  in  this 
application  to  verify  the  computed  velocities  or  the  implied  inundation 
extent.  It  would  be  desirable  to  consider  making  some  future  application 
of  the  model  for  reaches  where  such  data  are  available  to  explore  the 
interaction  that  may  exist  between  the  wetting/drying  parameters  used  in 
the  model  and  the  computed  overbank  velocities.  The  model  has  been  shown 
to  he  robust  against  field  uncertainty  in  floodplain  roughness 
estimation,  evidenced  by  the  results  shown  in  figures  7. 2-7. 5.  Stability 
of  solutions  for  wetting  and  drying  of  large  areas  was  greatly  improved 
by  use  or  the  "marsh"  element  option  in  version  4  of  RMA-2.  Further 
enhancement  of  this  capability,  as  well  as  that  of  channel  representation 
and  further  field  validation,  are  Identified  as  future  research  needs. 
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8.  CONCLUSIONS  AND  PROPOSALS  FOR  THE  NEXT  TWELVE  MONTHS 


8 . I  Conclusions 

1.  The  most  important  processes  active  in  two-stage  channels  flood 
prediction  have  been  identified  as  the  turbulent  exchange  between 
the  main  channel  and  floodplain  flow  segments  and  the 
short-circuiting  of  floodplain  flow  downstream. 

2.  These  two  processes  have  been  incorporated  in  MILHY3  whilst 
maintaining  parsimonious  data  requirements  and  retaining  the 
resolution  of  the  submodels. 

3.  For  events  where  at  least  15%  of  the  floodplain  is  apportioned  to 
the  floodplain,  the  application  of  the  multiple  routing  and 
turbulent  exchange  routines  together  improves  the  predictive 
capability  of  MILHY. 

4.  For  floodplain  where  less  than  15%  of  the  floodwave  is  apportioned 
to  the  floodplain,  only  the  multiple  routing  routine  has  any 
impact  on  the  hydrograph. 

5.  'Application  of  RMA2  shows  that  the  model  can  be  reliable  applied 

to  large  scale  reaches  and  that  it  could  provide  detailed 
inundation  predictions  for  ungauged  catchments,  utilizing  MILHY3 
to  provide  the  upstream  inflow  hydrograph. 

6.  Optimization  techniques  provide  a  viable  alternative  method  of 
tackling  factor  perturbation  sensitivity  analysis. 


8 . 2  Proposals  for  the  next  twelve  months 


1.  Continuing  investigation  of  the  feasibility  of  linking  MILHY3  with 
RMA-2  for  detailed  two-stage  channel  predictions  in  ungauged 
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catchments . 

Investigation  of  the  use  of  MILHY3/RMA-2  scheme  to  predict  the 
floodplain  erosion  and  sedimenattion  patterns.  This  investigation 
will  utilize  data  available  from  the  River  Culm,  Devon,  England. 
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